Load packages and data
options("lubridate.week.start" = 3)
library(tidyverse)
theme_set(theme_bw())
library(ggpubr)
library(lubridate)
library(ggrepel)
library(phyloseq)
library(grid)
library(ampvis2)
`%notin%` = Negate(`%in%`)
# Prevalence table function
# get overview of abundances, mean prevalence is the mean 'appearance' of ASVs of the taxon across all samples
prevalencedf <- dget("./custom-functions/myprevalencetablefunction.R")
pca_helper <- dget("./custom-functions/pca_helper.R")
startoperation <- ymd("2023-03-01")
cols <- c("#4D98AC", "#985C64", "#ff0000", "purple","#C4A484")
names(cols) <- c("Control", "Treatment", "Full-scale", "PSTWAS", "Foam")
## LOAD masterdata
# change this path to where you store the .csv file. you can provide an absolute path like this. Check the path syntax for Windows as it might be different.
masterdata <- read.csv("./data/masterdataProject1C_Jan25.csv")
masterdata <- masterdata %>%
mutate(Date = lubridate::dmy(as.character(Date)) )
# Make AD, Treatment and SludgeType = factors
masterdata$AD <- factor(masterdata$AD, levels = c("AD1", "AD2","AD3", "AD4","AD5", "AD6", "Full-scale", "PSTWAS"))
masterdata$SludgeType <- factor(masterdata$SludgeType, levels = c("Control", "Treatment", "Foam", "Full-scale", "PSTWAS"))
masterdata$Treatment <- factor(masterdata$Treatment, levels = c("Control", "Treatment", "Full-scale", "PSTWAS"))
masterdata$Period <- factor(masterdata$Period, levels = c("Converging", "SteadyState", "Glycerol", "Inhibition", "Recovery/Feedingpause", "Recovery/Foaming", "Recovery/Postfoam", "SteadyState/Postfoam"))
masterdata <- masterdata %>%
mutate(across(Observed:pielou_ev, ~as.numeric(.)))
# Phyloseq object with 16S abundance data
ps_1C <- readRDS("./data/ps_215samplesJul24") # Midas53
sample_data(ps_1C)$CSTR <- str_replace_all(sample_data(ps_1C)$AD, c("AD1"="C1",
"AD2"="T1",
"AD3"="C2",
"AD4"="C3",
"AD5"="T2",
"AD6"="T3",
"AD_full" = "Full-scale") )
ps.flt <- prune_taxa(taxa_sums(ps_1C) >= 10, ps_1C) #minimum reads per feature
prev.genus <- prevalencedf(ps.flt, Genus)
# Phyloseq object with metagenome abundance data
psmg <- readRDS("./data/ps-metagenome16samples")
## METADATA FOR METAGENOMICS SAMPLES
# merge individual samples into the pooled samples - mean values
df <- masterdata %>%
as_tibble() %>% dplyr::filter(SampleID.DNA %in% c("S1001_AD1", "S1003_AD3", "S1004_AD4",
"S1002_AD2","S1005_AD5","S1006_AD6",
"S1147_AD1", "S1149_AD3", "S1150_AD4",
"S1148_AD2", "S1151_AD5", "S1152_AD6",
"S1154_AD2","S1157_AD5","S1158_AD6"
)) %>% column_to_rownames("SampleID.DNA")
df <- data.frame(df %>% group_by(Date, SludgeType ) %>% summarise_each(funs(mean)))
rownames(df) <- c("S1001_3_4_AD1_3_4", "S1002_5_6_AD2_5_6","S1147_49_50_AD1_3_4", "S1148_51_52_AD2_5_6", "S1154_57_58_AD2_5_6")
df$Period <- c("Inhibition","Inhibition","Recovery/Foaming","Recovery/Foaming", "Recovery/Foaming")
df$AD <- c("AD1.3.4","AD2.5.6","AD1.3.4","AD2.5.6", "AD2.5.6")
# there was also 1 samples pooled for all 6 digesters
df1 <- masterdata %>%
as_tibble() %>% dplyr::filter(SampleID.DNA %in% c("S561_AD1", "S562_AD2", "S563_AD3",
"S564_AD4","S565_AD5","S566_AD6"
)) %>% column_to_rownames("SampleID.DNA")
df1 <- data.frame(df1 %>% group_by(Date ) %>% summarise_each(funs(mean)))
rownames(df1) <- c("S561toS566_AD.AD6")
df1$Period <- c("SteadyState")
df1$AD <- c("AD1.2.3.4.5.6")
# not pooled samples
df2 <- masterdata %>%
as_tibble() %>% dplyr::filter(SampleID.DNA %in% c("F5","F55","F110","F124","F179",
"S119_full","S370_full", "S1273_ADFull", "S1554_ADFull", "S1813_ADFull"
)) %>%
column_to_rownames("SampleID.DNA")
#combine the dfs
metadata <- rbind(df,df1, df2)
#rownames(metadata)
metadata <- metadata %>% rownames_to_column("ID") %>%
mutate("SampleID" = ID) %>% column_to_rownames("ID") %>%
mutate("SampleID" = factor(SampleID, levels = c("S561toS566_AD.AD6", "S1001_3_4_AD1_3_4", "S1002_5_6_AD2_5_6", "S1147_49_50_AD1_3_4", "S1148_51_52_AD2_5_6", "S1154_57_58_AD2_5_6", "F5", "F55", "F110","F124", "F179", "S119_full", "S370_full","S1273_ADFull", "S1554_ADFull","S1813_ADFull")))
metadata$ramaciottieID <- str_replace_all(metadata$SampleID, c( "F5"="Chris_01", "F55"="Chris_02", "F110"="Chris_03","F124" = "Chris_04",
"F179" = "Chris_05", "S119.full" = "Chris_06", "S370.full" = "Chris_07", "S1273_ADFull" = "Chris_08",
"S1554_ADFull" = "Chris_09", "S1813_ADFull" = "Chris_10", "S561toS566_AD.AD6" = "Chris_11", "S1001_3_4_AD1_3_4" = "Chris_12",
"S1002_5_6_AD2_5_6" = "Chris_13", "S1147_49_50_AD1_3_4" = "Chris_14", "S1148_51_52_AD2_5_6" = "Chris_15", "S1154_57_58_AD2_5_6" = "Chris_16"))
metadata <- metadata %>% mutate(VFAs = rowSums(dplyr::select(., contains("acid"))) )
metadata[c("S561toS566_AD.AD6"),c("SludgeType")] <- "Control"
rm(df, df1, df2)
## KEGG Pathways to PHYLOSEQ
## Load KO / KEGG pathway files ps object
ko_to_kegg_reference <- read.csv("./data/ko_to_kegg_reference.csv")[-1]
KO_reference <- read.csv("./data/KO_reference.csv")[-1]
KO_reference2 <- KO_reference %>% dplyr::select(-KO, -KoDescription) %>%
unique() %>%
separate_wider_delim(Pathway,
delim = ":",
names = c("Pathway","ko"),
too_few = "align_end") %>%
filter(Pathway != is.na(Pathway) ) %>%
mutate_at("ko", str_replace, "]", "") %>%
mutate_at("Pathway", str_replace, "PATH", "") %>%
mutate_at("Pathway", str_replace, "\\[", "")
KO_reference2$Pathway <- str_trim(KO_reference2$Pathway, "right")
KO_reference2 <- KO_reference2 %>% column_to_rownames("ko")
# abundances
KO_abundance <- read.csv("./data/KO.abundances_metagenome.csv")[-1]
KO_abundance <- KO_abundance %>% column_to_rownames("KO") %>%
dplyr::filter(rowSums(.) > 0) %>%
rownames_to_column("KO")
kegg_abundance <- read.csv("./data/kegg_abundances_metagenome.csv")[-1]
# Create the phyloseq object
ps_kegg <- phyloseq(
otu_table(kegg_abundance %>% column_to_rownames("ko"), taxa_are_rows = T),
sample_data(metadata),
phyloseq::tax_table(as.matrix(KO_reference2))
)
## PICRUST PATHWAYS TO PHYLOSEQ
# "otu_table"
pws <- data.frame(read_tsv("./data/path_abun_unstrat.tsv", col_names = TRUE)) %>%
rename(Id = "pathway")
# "tax_table"
MC_names <- read_tsv("./data/metacyc_pathway_info_full.tsv") %>%
dplyr::select(-ALTERNATIVE) %>%
dplyr::select(-`SUPER-PATHWAYS`)%>%
dplyr::select(-NAMES) %>%
dplyr::rename(Pathway2 = "COMMON-NAME")
# "tax_table new"
MC_names_new <- read_tsv("./data/metacyc_pathway_info_23jan.tsv")
#Combine
MC_names <- MC_names_new %>% dplyr::full_join(MC_names, by = "Id") %>%
filter(Id %in% pws$Id) %>%
column_to_rownames("Id")
# metadata
SplIDs <- rownames(sample_data(readRDS("./data/ps_215samplesJul24"))) %>% str_replace_all(".AD", "_AD")
SplIDs <- SplIDs %>% str_replace_all(".f", "_f")
metadata.pw <- masterdata %>%
as_tibble() %>% dplyr::filter(SampleID.DNA %in% SplIDs) %>% column_to_rownames("SampleID.DNA")
row.names(metadata.pw) <- row.names(metadata.pw) %>% str_replace_all("_AD", ".AD")
pws <- pws %>%
column_to_rownames("Id")
ps_pw <- phyloseq(
otu_table(pws, taxa_are_rows = T),
sample_data(metadata.pw),
phyloseq::tax_table(MC_names %>% as.matrix())
)
## KRAKEN to PHYLOSEQ
pskr <- readRDS("./data/pskr")
## ONP 16SFL to PHYLOSEQ
meta <- read.csv("./data/0824_sample_sheet16FL.csv") %>% column_to_rownames("alias")
meta$Sample_name <- factor(meta$Sample_name, levels = c( "Feed","Full-scale CSTR","CSTR1","CSTR2","CSTR3","CSTR4","CSTR5","CSTR6"))
meta$SampleType <- factor(meta$SampleType, levels = c( "Control","Treatment","Full-scale","PSTWAS"))
df <- read_tsv("./data/0824_abundance_table_species_16FL.tsv")
taxa <- df %>%
separate(tax, sep = ";", into = c("Superkingdom","Kingdom","Phylum","Class","Order","Family","Genus","Species"))%>%
dplyr::select(Superkingdom,Phylum,Kingdom,Class,Order,Family,Genus,Species) %>%
column_to_rownames("Species")
taxa$Species <- rownames(taxa)
otus <- df[,c(2:9)]
rownames(otus) <- taxa$Species
# Create the phyloseq object
ps.onp.ncbi.kraken <- phyloseq(
otu_table(otus, taxa_are_rows = T),
sample_data(meta),
phyloseq::tax_table(as.matrix(taxa))
)
## FOAM DATA
foamdata <- read.csv("./data/foam_potential_list3Copy_CK.csv")
foamdata <- foamdata %>%
mutate(Date = lubridate::dmy(as.character(Date)) )
foamdata <- foamdata %>% column_to_rownames("X.SludgeID")
Supp. figure S3 gflow
# Daily from masterfile
nrow(masterdata %>%
dplyr::filter(Gasflow_mL > 20000)) # outlier
## [1] 2
generate_vector <- function(x, y) {
result <- c()
while (x <= y) {
result <- c(result, x)
x <- x + 7
}
return(result)
}
generate_vectorday <- function(x, y) {
result <- c()
while (x <= y) {
result <- c(result, x)
x <- x + 1
}
return(result)
}
firstdate <- -2 # selected periods for Supplementary
lastdate <- 166
weeky <- -1
pexcl <- masterdata %>%
dplyr::filter(Gasflow_mL >= 0 &
Gasflow_mL < 20000) %>%
dplyr::filter(Date != c("2023-06-08")) %>% # outlier as the heaters turned off that day
dplyr::filter(AD != "AD4") %>% #AD4 always behaved differently due to the lack of pre-mixing
mutate(Gasflow_L = Gasflow_mL / 1000) %>%
ggline(x = "daysop",
y = "Gasflow_L",
color = "Treatment",
legend = "right",
add = "mean_se",
add.params = list(width = 0.75)) +
scale_color_manual("CSTR", values = cols) +
theme_light() +
theme(
#axis.title.x = element_blank(),
#axis.text.x = element_text(angle = 0, hjust=1),
legend.position = "none",
legend.background=element_rect(fill = alpha("white", 0.5)) ) +
scale_x_continuous(breaks = generate_vector(firstdate,lastdate),
minor_breaks = generate_vectorday(firstdate,lastdate),
limits = c(firstdate, lastdate) ) +
# scale_x_date(date_breaks = "1 week", date_labels = "%b %d",
# date_minor_breaks = "1 day",
# limits = as.Date(c(firstdate, lastdate)) ) +
annotate("text",x = 2, y = weeky, label = "W01 ", size = 2.5) + #Week 01
annotate("text",x = 9, y = weeky, label = "W02 ", size = 2.5) + #Week 02
annotate("text",x = 16, y = weeky, label = "W03 ", size = 2.5) + #Week 03
annotate("text",x = 23, y = weeky, label = "W04 ", size = 2.5) + #Week 04
annotate("text",x = 30, y = weeky, label = "W05 ", size = 2.5) + #Week 05
annotate("text",x = 37, y = weeky, label = "W06 ", size = 2.5) + #Week 06
annotate("text",x = 44, y = weeky, label = "W07 ", size = 2.5) + #Week 07
annotate("text",x = 51, y = weeky, label = "W08 ", size = 2.5) + #Week 08
annotate("text",x = 58, y = weeky, label = "W09 ", size = 2.5) + #Week 09
annotate("text",x = 65, y = weeky, label = "W10 ", size = 2.5) + #Week 10
annotate("text",x = 72, y = weeky, label = "W11 ", size = 2.5) + #Week 11
annotate("text",x = 79, y = weeky, label = "W12 ", size = 2.5) + #Week 12
annotate("text",x = 86, y = weeky, label = "W13 ", size = 2.5) + #Week 13
annotate("text",x = 93, y = weeky, label = "W14 ", size = 2.5) + #Week 14
annotate("text",x = 100, y = weeky, label = "W15 ", size = 2.5) + #Week 15
annotate("text",x = 107, y = weeky, label = "W16 ", size = 2.5) + #Week 16
annotate("text",x = 114, y = weeky, label = "W17 ", size = 2.5) + #Week 17
annotate("text",x = 121, y = weeky, label = "W18 ", size = 2.5) + #Week 18
annotate("text",x = 128, y = weeky, label = "W19 ", size = 2.5) + #Week 19
annotate("text",x = 135, y = weeky, label = "W20 ", size = 2.5) + #Week 20
annotate("text",x = 142, y = weeky, label = "W21 ", size = 2.5) + #Week 21
annotate("text",x = 149, y = weeky, label = "W22 ", size = 2.5) + #Week 22
annotate("text",x = 156, y = weeky, label = "W23 ", size = 2.5) + #Week 23
annotate("text",x = 163, y = weeky, label = "W24 ", size = 2.5) + #Week 24
annotate("text",x = 19, y = -.4, label = "Converging", size = 2.5) +
annotate("segment", y = 0, yend = 0, x = 0, xend = 31,
arrow = arrow(ends = "last", angle = 90, length = unit(.1,"cm"))) +
annotate("text",x = 47, y = -.4, label = "Steady state", size = 2.5) +
annotate("segment", y = 0, yend = 0, x = 31, xend = 62,
arrow = arrow(ends = "last", angle = 90, length = unit(.1,"cm"))) +
annotate("text",x = 65, y = -.4, label = "Glycerol ", size = 2.5) +
annotate("segment", y = 0, yend = 0, x = 62, xend = 67,
arrow = arrow(ends = "both", angle = 90, length = unit(.1,"cm"))) +
annotate("text",x = 74, y = -0.4, label = " Inhibition", size = 2.5) +
annotate("segment", y = 0, yend = 0, x = 67, xend = 82,
arrow = arrow(ends = "both", angle = 90, length = unit(.1,"cm"))) +
annotate("text",x = 86, y = -0.4, label = "Recovery", size = 2.5) +
annotate("segment", y = 0, yend = 0, x = 82, xend = 90,
arrow = arrow(ends = "both", angle = 90, length = unit(.1,"cm"))) +
annotate("text",x = 94, y = -0.4, label = "Foaming ", size = 2.5) +
annotate("segment", y = 0, yend = 0, x = 90, xend = 98,
arrow = arrow(ends = "both", angle = 90, length = unit(.1,"cm"))) +
annotate("text",x = 101, y = -0.4, label = " Recovery", size = 2.5) +
annotate("segment", y = 0, yend = 0, x = 98, xend = 106,
arrow = arrow(ends = "both", angle = 90, length = unit(.1,"cm"))) +
annotate("text",x = 138, y = -.4, label = "Steady state post foam", size = 3) +
annotate("segment", y = 0, yend = 0, x = 98, xend = 164,
arrow = arrow(ends = "first", angle = 90, length = unit(.1,"cm"))) +
annotate("text",x = 85, y = 3.65, label = "Feeding paused", size = 3, angle = 90, alpha = 0.6) +
annotate("rect", ymin = 0, ymax = 6, xmin = 82, xmax = 89,
fill = "blue", alpha = .1) +
annotate("text",x = 94, y = 7, label = "Foam", size = 3, angle = 90, alpha = 0.6) +
annotate("rect", ymin = 4.5, ymax = 8, xmin = 92, xmax = 97,
fill = "blue", alpha = .1) +
ylab(expression(paste("Gasflow (L day"^-1*")"))) +
xlab("Day of operation")
pincl <- masterdata %>%
dplyr::filter(Gasflow_mL >= 0 &
Gasflow_mL < 20000) %>%
dplyr::filter(Date != c("2023-06-08")) %>% # outlier as the heaters turned off that day
# dplyr::filter(AD != "AD4") %>% #AD4 always behaved differently due to the lack of pre-mixing
mutate(Gasflow_L = Gasflow_mL / 1000) %>%
ggline(x = "daysop",
y = "Gasflow_L",
color = "Treatment",
legend = "right",
add = "mean_se",
add.params = list(width = 0.75)) +
scale_color_manual("CSTR", values = cols) +
theme_light() +
theme(
axis.title.x = element_blank(),
#axis.text.x = element_text(angle = 0, hjust=1),
legend.position = c(0.90, 0.82),
legend.background=element_rect(fill = alpha("white", 0.5)) ) +
scale_x_continuous(breaks = generate_vector(firstdate,lastdate),
minor_breaks = generate_vectorday(firstdate,lastdate),
limits = c(firstdate, lastdate) ) +
# scale_x_date(date_breaks = "1 week", date_labels = "%b %d",
# date_minor_breaks = "1 day",
# limits = as.Date(c(firstdate, lastdate)) ) +
annotate("text",x = 2, y = weeky, label = "W01 ", size = 2.5) + #Week 01
annotate("text",x = 9, y = weeky, label = "W02 ", size = 2.5) + #Week 02
annotate("text",x = 16, y = weeky, label = "W03 ", size = 2.5) + #Week 03
annotate("text",x = 23, y = weeky, label = "W04 ", size = 2.5) + #Week 04
annotate("text",x = 30, y = weeky, label = "W05 ", size = 2.5) + #Week 05
annotate("text",x = 37, y = weeky, label = "W06 ", size = 2.5) + #Week 06
annotate("text",x = 44, y = weeky, label = "W07 ", size = 2.5) + #Week 07
annotate("text",x = 51, y = weeky, label = "W08 ", size = 2.5) + #Week 08
annotate("text",x = 58, y = weeky, label = "W09 ", size = 2.5) + #Week 09
annotate("text",x = 65, y = weeky, label = "W10 ", size = 2.5) + #Week 10
annotate("text",x = 72, y = weeky, label = "W11 ", size = 2.5) + #Week 11
annotate("text",x = 79, y = weeky, label = "W12 ", size = 2.5) + #Week 12
annotate("text",x = 86, y = weeky, label = "W13 ", size = 2.5) + #Week 13
annotate("text",x = 93, y = weeky, label = "W14 ", size = 2.5) + #Week 14
annotate("text",x = 100, y = weeky, label = "W15 ", size = 2.5) + #Week 15
annotate("text",x = 107, y = weeky, label = "W16 ", size = 2.5) + #Week 16
annotate("text",x = 114, y = weeky, label = "W17 ", size = 2.5) + #Week 17
annotate("text",x = 121, y = weeky, label = "W18 ", size = 2.5) + #Week 18
annotate("text",x = 128, y = weeky, label = "W19 ", size = 2.5) + #Week 19
annotate("text",x = 135, y = weeky, label = "W20 ", size = 2.5) + #Week 20
annotate("text",x = 142, y = weeky, label = "W21 ", size = 2.5) + #Week 21
annotate("text",x = 149, y = weeky, label = "W22 ", size = 2.5) + #Week 22
annotate("text",x = 156, y = weeky, label = "W23 ", size = 2.5) + #Week 23
annotate("text",x = 163, y = weeky, label = "W24 ", size = 2.5) + #Week 24
annotate("text",x = 19, y = -.4, label = "Converging", size = 2.5) +
annotate("segment", y = 0, yend = 0, x = 0, xend = 31,
arrow = arrow(ends = "last", angle = 90, length = unit(.1,"cm"))) +
annotate("text",x = 47, y = -.4, label = "Steady state", size = 2.5) +
annotate("segment", y = 0, yend = 0, x = 31, xend = 62,
arrow = arrow(ends = "last", angle = 90, length = unit(.1,"cm"))) +
annotate("text",x = 65, y = -.4, label = "Glycerol ", size = 2.5) +
annotate("segment", y = 0, yend = 0, x = 62, xend = 67,
arrow = arrow(ends = "both", angle = 90, length = unit(.1,"cm"))) +
annotate("text",x = 74, y = -0.4, label = " Inhibition", size = 2.5) +
annotate("segment", y = 0, yend = 0, x = 67, xend = 82,
arrow = arrow(ends = "both", angle = 90, length = unit(.1,"cm"))) +
annotate("text",x = 86, y = -0.4, label = "Recovery", size = 2.5) +
annotate("segment", y = 0, yend = 0, x = 82, xend = 90,
arrow = arrow(ends = "both", angle = 90, length = unit(.1,"cm"))) +
annotate("text",x = 94, y = -0.4, label = "Foaming ", size = 2.5) +
annotate("segment", y = 0, yend = 0, x = 90, xend = 98,
arrow = arrow(ends = "both", angle = 90, length = unit(.1,"cm"))) +
annotate("text",x = 101, y = -0.4, label = " Recovery", size = 2.5) +
annotate("segment", y = 0, yend = 0, x = 98, xend = 106,
arrow = arrow(ends = "both", angle = 90, length = unit(.1,"cm"))) +
annotate("text",x = 138, y = -.4, label = "Steady state post foam", size = 3) +
annotate("segment", y = 0, yend = 0, x = 98, xend = 164,
arrow = arrow(ends = "first", angle = 90, length = unit(.1,"cm"))) +
annotate("text",x = 85, y = 3.65, label = "Feeding paused", size = 3, angle = 90, alpha = 0.6) +
annotate("rect", ymin = 0, ymax = 6, xmin = 82, xmax = 89,
fill = "blue", alpha = .1) +
annotate("text",x = 94, y = 7, label = "Foam", size = 3, angle = 90, alpha = 0.6) +
annotate("rect", ymin = 4.5, ymax = 8, xmin = 92, xmax = 97,
fill = "blue", alpha = .1) +
ylab(expression(paste("Gasflow (L day"^-1*")")))
ggarrange(pincl, pexcl, labels = "AUTO", ncol = 1)
ggsave("./Figures/Sup_lineplot_gasdaily_inclexclAD4.png", height=20, width=27, units='cm', dpi=300)
Supp. figure S3 gflow - ind ADs (not used)
# Daily from masterfile
nrow(masterdata %>%
dplyr::filter(Gasflow_mL > 20000)) # outlier
## [1] 2
firstdate <- "2023-03-01"
lastdate <- "2023-08-15"
masterdata %>%
dplyr::filter(Gasflow_mL >= 0 &
Gasflow_mL < 20000) %>%
dplyr::filter(Date != c("2023-06-08")) %>% # outlier as the heaters turned off that day
# dplyr::filter(AD != "AD4") %>% #AD4 always behaved differently due to the lack of pre-mixing
mutate(Gasflow_L = Gasflow_mL / 1000) %>%
ggline(x = "Date",
y = "Gasflow_L",
color = "AD",
legend = "right",
add = "mean_se",
add.params = list(width = 0.75)) +
#scale_color_manual("Sludge Type", values = cols) +
theme_light() +
theme(axis.title.x = element_blank(),
axis.text.x = element_text(angle = 45, hjust=1),
legend.position = c(0.90, 0.82),
legend.background=element_rect(fill = alpha("white", 0.5)) ) +
scale_x_date(date_breaks = "1 week", date_labels = "%b %d",
date_minor_breaks = "1 day",
limits = as.Date(c(firstdate, lastdate)) ) +
annotate("text",x = as.Date("2023-04-21"), y = -1, label = "W08 ", size = 2.5) + #Week 08
annotate("text",x = as.Date("2023-04-28"), y = -1, label = "W09 ", size = 2.5) + #Week 09
annotate("text",x = as.Date("2023-05-05"), y = -1, label = "W10 ", size = 2.5) + #Week 10
annotate("text",x = as.Date("2023-05-12"), y = -1, label = "W11 ", size = 2.5) + #Week 11
annotate("text",x = as.Date("2023-05-19"), y = -1, label = "W12 ", size = 2.5) + #Week 12
annotate("text",x = as.Date("2023-05-26"), y = -1, label = "W13 ", size = 2.5) + #Week 13
annotate("text",x = as.Date("2023-06-02"), y = -1, label = "W14 ", size = 2.5) + #Week 14
annotate("text",x = as.Date("2023-06-09"), y = -1, label = "W15 ", size = 2.5) + #Week 15
annotate("text",x = as.Date("2023-06-16"), y = -1, label = "W16 ", size = 2.5) + #Week 16
annotate("text",x = as.Date("2023-06-23"), y = -1, label = "W17 ", size = 2.5) + #Week 17
annotate("text",x = as.Date("2023-06-30"), y = -1, label = "W18 ", size = 2.5) + #Week 18
annotate("text",x = as.Date("2023-04-24"), y = -.4, label = "Steady state", size = 2.5) +
annotate("segment", y = 0, yend = 0, x = as.Date("2023-04-17"), xend = as.Date("2023-05-02"),
arrow = arrow(ends = "last", angle = 90, length = unit(.1,"cm"))) +
annotate("text",x = as.Date("2023-05-05"), y = -.4, label = "Glycerol ", size = 2.5) +
annotate("segment", y = 0, yend = 0, x = as.Date("2023-05-02"), xend = as.Date("2023-05-07"),
arrow = arrow(ends = "both", angle = 90, length = unit(.1,"cm"))) +
annotate("text",x = as.Date("2023-05-14"), y = -0.4, label = " Inhibition", size = 2.5) +
annotate("segment", y = 0, yend = 0, x = as.Date("2023-05-07"), xend = as.Date("2023-05-22"),
arrow = arrow(ends = "both", angle = 90, length = unit(.1,"cm"))) +
annotate("text",x = as.Date("2023-05-26"), y = -0.4, label = "Recovery", size = 2.5) +
annotate("segment", y = 0, yend = 0, x = as.Date("2023-05-22"), xend = as.Date("2023-05-30"),
arrow = arrow(ends = "both", angle = 90, length = unit(.1,"cm"))) +
annotate("text",x = as.Date("2023-06-03"), y = -0.4, label = "Foaming ", size = 2.5) +
annotate("segment", y = 0, yend = 0, x = as.Date("2023-05-30"), xend = as.Date("2023-06-07"),
arrow = arrow(ends = "both", angle = 90, length = unit(.1,"cm"))) +
annotate("text",x = as.Date("2023-06-10"), y = -0.4, label = " Recovery", size = 2.5) +
annotate("segment", y = 0, yend = 0, x = as.Date("2023-06-07"), xend = as.Date("2023-06-15"),
arrow = arrow(ends = "both", angle = 90, length = unit(.1,"cm"))) +
annotate("text",x = as.Date("2023-05-25"), y = 4, label = "Paused feeding", size = 2.5, angle = 90, alpha = 0.6) +
annotate("rect", ymin = 0, ymax = 6, xmin = as.Date("2023-05-22"), xmax = as.Date("2023-05-29"),
fill = "blue", alpha = .1) +
annotate("text",x = as.Date("2023-06-03"), y = 7, label = "Foam", size = 2.5, angle = 90, alpha = 0.6) +
annotate("rect", ymin = 4.5, ymax = 8, xmin = as.Date("2023-06-01"), xmax = as.Date("2023-06-06"),
fill = "blue", alpha = .1) +
annotate("text",x = as.Date("2023-06-24"), y = -.4, label = "Steady state post foam", size = 2.5) +
annotate("segment", y = 0, yend = 0, x = as.Date("2023-06-07"), xend = as.Date("2023-07-03"),
arrow = arrow(ends = "first", angle = 90, length = unit(.1,"cm"))) +
ylab(expression(paste("Gasflow (L day"^-1*")")))
ggsave("./Figures/Sup_lineplot_gasdaily_indADs.png", height=15, width=25, units='cm', dpi=300)
Supp. figure S4 CH4
firstdate <- -2 # selected periods for Supplementary
lastdate <- 166
weeky <- c(-6.5)
labely <- c(-2.5)
segmenty <- c(0)
generate_vector <- function(x, y) {
result <- c()
while (x <= y) {
result <- c(result, x)
x <- x + 7
}
return(result)
}
generate_vectorday <- function(x, y) {
result <- c()
while (x <= y) {
result <- c(result, x)
x <- x + 1
}
return(result)
}
## Gas comp
masterdata %>%
dplyr::filter(CH4_pct >= 0) %>%
ggline(x = "daysop",
y = "CH4_pct",
color = "Treatment",
legend = "right",
add = "mean_se",
add.params = list(width = 1)) +
ylab("CH4 (%)") +
theme_light() +
scale_color_manual("CSTR", values = cols) +
theme(
#axis.title.x = element_blank(),
# axis.text.x = element_text(angle = 0, hjust=1),
legend.position = c(0.9, 0.4)
# legend.position = "none"
) +
scale_x_continuous(breaks = generate_vector(firstdate,lastdate),
minor_breaks = generate_vectorday(firstdate,lastdate),
limits = c(firstdate, lastdate) ) +
annotate("text",x = 2, y = weeky, label = "W01 ", size = 2.5) + #Week 01
annotate("text",x = 9, y = weeky, label = "W02 ", size = 2.5) + #Week 02
annotate("text",x = 16, y = weeky, label = "W03 ", size = 2.5) + #Week 03
annotate("text",x = 23, y = weeky, label = "W04 ", size = 2.5) + #Week 04
annotate("text",x = 30, y = weeky, label = "W05 ", size = 2.5) + #Week 05
annotate("text",x = 37, y = weeky, label = "W06 ", size = 2.5) + #Week 06
annotate("text",x = 44, y = weeky, label = "W07 ", size = 2.5) + #Week 07
annotate("text",x = 51, y = weeky, label = "W08 ", size = 2.5) + #Week 08
annotate("text",x = 58, y = weeky, label = "W09 ", size = 2.5) + #Week 09
annotate("text",x = 65, y = weeky, label = "W10 ", size = 2.5) + #Week 10
annotate("text",x = 72, y = weeky, label = "W11 ", size = 2.5) + #Week 11
annotate("text",x = 79, y = weeky, label = "W12 ", size = 2.5) + #Week 12
annotate("text",x = 86, y = weeky, label = "W13 ", size = 2.5) + #Week 13
annotate("text",x = 93, y = weeky, label = "W14 ", size = 2.5) + #Week 14
annotate("text",x = 100, y = weeky, label = "W15 ", size = 2.5) + #Week 15
annotate("text",x = 107, y = weeky, label = "W16 ", size = 2.5) + #Week 16
annotate("text",x = 114, y = weeky, label = "W17 ", size = 2.5) + #Week 17
annotate("text",x = 121, y = weeky, label = "W18 ", size = 2.5) + #Week 18
annotate("text",x = 128, y = weeky, label = "W19 ", size = 2.5) + #Week 19
annotate("text",x = 135, y = weeky, label = "W20 ", size = 2.5) + #Week 20
annotate("text",x = 142, y = weeky, label = "W21 ", size = 2.5) + #Week 21
annotate("text",x = 149, y = weeky, label = "W22 ", size = 2.5) + #Week 22
annotate("text",x = 156, y = weeky, label = "W23 ", size = 2.5) + #Week 23
annotate("text",x = 163, y = weeky, label = "W24 ", size = 2.5) + #Week 24
annotate("text",x = 19, y = labely, label = "Converging", size = 2.5) +
annotate("segment", y = 0, yend = 0, x = 0, xend = 31,
arrow = arrow(ends = "last", angle = 90, length = unit(.1,"cm"))) +
annotate("text",x = 47, y = labely, label = "Steady state", size = 2.5) +
annotate("segment", y = 0, yend = 0, x = 31, xend = 62,
arrow = arrow(ends = "last", angle = 90, length = unit(.1,"cm"))) +
annotate("text",x = 65, y = labely, label = "Glycerol ", size = 2.5) +
annotate("segment", y = 0, yend = 0, x = 62, xend = 67,
arrow = arrow(ends = "both", angle = 90, length = unit(.1,"cm"))) +
annotate("text",x = 74, y = labely, label = " Inhibition", size = 2.5) +
annotate("segment", y = 0, yend = 0, x = 67, xend = 82,
arrow = arrow(ends = "both", angle = 90, length = unit(.1,"cm"))) +
annotate("text",x = 86, y = labely, label = "Recovery", size = 2.5) +
annotate("segment", y = 0, yend = 0, x = 82, xend = 90,
arrow = arrow(ends = "both", angle = 90, length = unit(.1,"cm"))) +
annotate("text",x = 94, y = labely, label = "Foaming ", size = 2.5) +
annotate("segment", y = 0, yend = 0, x = 90, xend = 98,
arrow = arrow(ends = "both", angle = 90, length = unit(.1,"cm"))) +
annotate("text",x = 101, y = labely, label = " Recovery", size = 2.5) +
annotate("segment", y = 0, yend = 0, x = 98, xend = 106,
arrow = arrow(ends = "both", angle = 90, length = unit(.1,"cm"))) +
annotate("text",x = 138, y = labely, label = "Steady state post foam", size = 2.5) +
annotate("segment", y = 0, yend = 0, x = 98, xend = 164,
arrow = arrow(ends = "first", angle = 90, length = unit(.1,"cm"))) +
annotate("text",x = 85.5, y = 35, label = "Feeding\npaused", size = 3.5, alpha = 0.4, angle = 90) +
annotate("rect", ymin = 0, ymax = 75, xmin = 82, xmax = 89,
fill = "blue", alpha = .1) +
annotate("text",x = 94.5, y = 67, label = "Foaming", size = 2.5, alpha = 0.5, angle = 90) +
annotate("rect", ymin = 60, ymax = 75, xmin = 92, xmax = 97, fill = "blue", alpha = .1) +
#
ylab(expression(paste("CH"[4]*" (%)"))) +
xlab("Day of operation")
ggsave("./Figures/Sup_lineplot_CH4.png", height=8, width=26, units='cm', dpi=300)
Supp. figure S5 pH
firstdate <- -2 # selected periods for Supplementary
lastdate <- 166
weeky <- c(4.8)
labely <- c(5.25)
segmenty <- c(5)
generate_vector <- function(x, y) {
result <- c()
while (x <= y) {
result <- c(result, x)
x <- x + 7
}
return(result)
}
generate_vectorday <- function(x, y) {
result <- c()
while (x <= y) {
result <- c(result, x)
x <- x + 1
}
return(result)
}
masterdata %>%
dplyr::filter(Treatment != "PSTWAS") %>%
dplyr::filter(pH >= 0) %>%
ggline(x = "daysop",
y = "pH",
color = "Treatment",
legend = "right",
add = "mean_se",
add.params = list(width = 1)) +
scale_color_manual("CSTR", values = cols) +
theme_light() +
theme(
#axis.title.x = element_blank(),
#axis.text.x = element_text(angle = 0, hjust=1),
legend.position = c(0.9, 0.4)) +
scale_x_continuous(breaks = generate_vector(firstdate,lastdate),
minor_breaks = generate_vectorday(firstdate,lastdate),
limits = c(firstdate, lastdate) ) +
annotate("text",x = 2, y = weeky, label = "W01 ", size = 2.5) + #Week 01
annotate("text",x = 9, y = weeky, label = "W02 ", size = 2.5) + #Week 02
annotate("text",x = 16, y = weeky, label = "W03 ", size = 2.5) + #Week 03
annotate("text",x = 23, y = weeky, label = "W04 ", size = 2.5) + #Week 04
annotate("text",x = 30, y = weeky, label = "W05 ", size = 2.5) + #Week 05
annotate("text",x = 37, y = weeky, label = "W06 ", size = 2.5) + #Week 06
annotate("text",x = 44, y = weeky, label = "W07 ", size = 2.5) + #Week 07
annotate("text",x = 51, y = weeky, label = "W08 ", size = 2.5) + #Week 08
annotate("text",x = 58, y = weeky, label = "W09 ", size = 2.5) + #Week 09
annotate("text",x = 65, y = weeky, label = "W10 ", size = 2.5) + #Week 10
annotate("text",x = 72, y = weeky, label = "W11 ", size = 2.5) + #Week 11
annotate("text",x = 79, y = weeky, label = "W12 ", size = 2.5) + #Week 12
annotate("text",x = 86, y = weeky, label = "W13 ", size = 2.5) + #Week 13
annotate("text",x = 93, y = weeky, label = "W14 ", size = 2.5) + #Week 14
annotate("text",x = 100, y = weeky, label = "W15 ", size = 2.5) + #Week 15
annotate("text",x = 107, y = weeky, label = "W16 ", size = 2.5) + #Week 16
annotate("text",x = 114, y = weeky, label = "W17 ", size = 2.5) + #Week 17
annotate("text",x = 121, y = weeky, label = "W18 ", size = 2.5) + #Week 18
annotate("text",x = 128, y = weeky, label = "W19 ", size = 2.5) + #Week 19
annotate("text",x = 135, y = weeky, label = "W20 ", size = 2.5) + #Week 20
annotate("text",x = 142, y = weeky, label = "W21 ", size = 2.5) + #Week 21
annotate("text",x = 149, y = weeky, label = "W22 ", size = 2.5) + #Week 22
annotate("text",x = 156, y = weeky, label = "W23 ", size = 2.5) + #Week 23
annotate("text",x = 163, y = weeky, label = "W24 ", size = 2.5) + #Week 24
annotate("text",x = 19, y = labely, label = "Converging", size = 2.5) +
annotate("segment", y = segmenty, yend = segmenty, x = 0, xend = 31,
arrow = arrow(ends = "last", angle = 90, length = unit(.1,"cm"))) +
annotate("text",x = 47, y = labely, label = "Steady state", size = 2.5) +
annotate("segment", y = segmenty, yend = segmenty, x = 31, xend = 62,
arrow = arrow(ends = "last", angle = 90, length = unit(.1,"cm"))) +
annotate("text",x = 65, y = labely, label = "Glycerol ", size = 2.5) +
annotate("segment", y = segmenty, yend = segmenty, x = 62, xend = 67,
arrow = arrow(ends = "both", angle = 90, length = unit(.1,"cm"))) +
annotate("text",x = 74, y = labely, label = " Inhibition", size = 2.5) +
annotate("segment", y = segmenty, yend = segmenty, x = 67, xend = 82,
arrow = arrow(ends = "both", angle = 90, length = unit(.1,"cm"))) +
annotate("text",x = 86, y = labely, label = "Recovery", size = 2.5) +
annotate("segment", y = segmenty, yend = segmenty, x = 82, xend = 90,
arrow = arrow(ends = "both", angle = 90, length = unit(.1,"cm"))) +
annotate("text",x = 94, y = labely, label = "Foaming ", size = 2.5) +
annotate("segment", y = segmenty, yend = segmenty, x = 90, xend = 98,
arrow = arrow(ends = "both", angle = 90, length = unit(.1,"cm"))) +
annotate("text",x = 101, y = labely, label = " Recovery", size = 2.5) +
annotate("segment", y = segmenty, yend = segmenty, x = 98, xend = 106,
arrow = arrow(ends = "both", angle = 90, length = unit(.1,"cm"))) +
annotate("text",x = 138, y = labely, label = "Steady state post foam", size = 2.5) +
annotate("segment", y = segmenty, yend = segmenty, x = 98, xend = 164,
arrow = arrow(ends = "first", angle = 90, length = unit(.1,"cm"))) +
annotate("text",x = 85.5, y = 7, label = "Feeding\npaused", size = 3.5, alpha = 0.4, angle = 90) +
annotate("rect", ymin = 6, ymax = 7.75, xmin = 82, xmax = 89,
fill = "blue", alpha = .1) +
ylab("pH") +
xlab("Day of operation")
ggsave("./Figures/Sup_lineplot_ph.png", height=8, width=26, units='cm', dpi=300)
Supp. figure S6 DNA (mg gCOD)
firstdate <- 47
lastdate <- 124
weeky <- -0.3
labely <- c(-0.125)
segmenty <- c(0)
masterdata %>% dplyr::filter(DNA_mg_gCOD > 0 &
DNA_mg_gCOD < 5 ) %>%
ggline(x = "daysop",
y = "DNA_mg_gCOD",
color = "SludgeType",
legend = "right",
add = "mean_se",
add.params = list(width = 1)) +
scale_color_manual("Sludge type",
values = cols, labels=c('Control CSTR', 'Treatment CSTR', 'Full-scale CSTR', "PS/TWAS")) +
scale_shape_manual(values = c(15, 15, 16, 17)) +
guides(color = guide_legend(
override.aes=list(shape = c(15, 15, 16, 17))),
shape = "none") +
theme_light() +
theme(
#axis.title.x = element_blank(),
#axis.text.x = element_text(angle = 0, hjust=.5)
legend.position = c(0.88, 0.83),
legend.background=element_rect(fill = alpha("white", 0.1)),
legend.key.height = unit(0.7, "line") # Adjust this value to make the legend shorter
) +
scale_x_continuous(breaks = generate_vector(firstdate,lastdate),
minor_breaks = generate_vectorday(firstdate,lastdate),
limits = c(firstdate, lastdate) ) +
annotate("text",x = 51, y = weeky, label = "W08 ", size = 2.5) + #Week 08
annotate("text",x = 58, y = weeky, label = "W09 ", size = 2.5) + #Week 09
annotate("text",x = 65, y = weeky, label = "W10 ", size = 2.5) + #Week 10
annotate("text",x = 72, y = weeky, label = "W11 ", size = 2.5) + #Week 11
annotate("text",x = 79, y = weeky, label = "W12 ", size = 2.5) + #Week 12
annotate("text",x = 86, y = weeky, label = "W13 ", size = 2.5) + #Week 13
annotate("text",x = 93, y = weeky, label = "W14 ", size = 2.5) + #Week 14
annotate("text",x = 100, y = weeky, label = "W15 ", size = 2.5) + #Week 15
annotate("text",x = 107, y = weeky, label = "W16 ", size = 2.5) + #Week 16
annotate("text",x = 114, y = weeky, label = "W17 ", size = 2.5) + #Week 17
annotate("text",x = 121, y = weeky, label = "W18 ", size = 2.5) + #Week 18
annotate("text",x = 54, y = labely, label = "Steady state", size = 2.5) +
annotate("segment", y = segmenty, yend = segmenty, x = 47, xend = 62,
arrow = arrow(ends = "last", angle = 90, length = unit(.1,"cm"))) +
annotate("text",x = 65, y = labely, label = "Glycerol ", size = 2.5) +
annotate("segment", y = segmenty, yend = segmenty, x = 62, xend = 67,
arrow = arrow(ends = "both", angle = 90, length = unit(.1,"cm"))) +
annotate("text",x = 74, y = labely, label = " Inhibition", size = 2.5) +
annotate("segment", y = segmenty, yend = segmenty, x = 67, xend = 82,
arrow = arrow(ends = "both", angle = 90, length = unit(.1,"cm"))) +
annotate("text",x = 85, y = labely, label = " Recovery", size = 2.5) +
annotate("segment", y = segmenty, yend = segmenty, x = 82, xend = 90,
arrow = arrow(ends = "both", angle = 90, length = unit(.1,"cm"))) +
annotate("text",x = 94, y = labely, label = "Foaming ", size = 2.5) +
annotate("segment", y = segmenty, yend = segmenty, x = 90, xend = 98,
arrow = arrow(ends = "both", angle = 90, length = unit(.1,"cm"))) +
annotate("text",x = 101, y = labely, label = " Recovery", size = 2.5) +
annotate("segment", y = segmenty, yend = segmenty, x = 98, xend = 106,
arrow = arrow(ends = "both", angle = 90, length = unit(.1,"cm"))) +
annotate("text",x = 115, y = labely, label = "Steady state post foam", size = 2.5) +
annotate("segment", y = segmenty, yend = segmenty, x = 98, xend = 124,
arrow = arrow(ends = "first", angle = 90, length = unit(.1,"cm"))) +
annotate("text",x = 85, y = 1.25, label = "Feeding paused", size = 3, angle = 90, alpha = 0.5) +
annotate("rect", ymin = 0, ymax = 2.5, xmin = 82, xmax = 89, fill = "blue", alpha = .1) +
ylab(expression(paste("DNA (mg gCOD"^-1*")"))) +
xlab("Day of operation")
ggsave("./Figures/Sup_lineplot_mgCOD.png", height=3, width=6.75, units='in', dpi=300)
Supp. figure S7 Shannon
firstdate <- 47
lastdate <- 124
weeky <- 3.2
labely <- c(3.35)
segmenty <- c(3.5)
masterdata %>% dplyr::filter(ngmL_qb >= 0 ) %>%
ggline(x = "daysop",
y = "Shannon",
color = "SludgeType",
legend = "top",
add = "mean_se",
add.params = list(width = 1)) +
scale_color_manual("Sample type",
values = cols, labels=c('Control CSTR', 'Treatment CSTR', 'Foam', "Full-scale CSTR", "PS/TWAS")) +
scale_shape_manual(values = c(15, 15, 16, 15, 17)) +
guides(color = guide_legend(
override.aes=list(shape = c(15, 15, 16, 15, 17))),
shape = "none") +
theme_light() +
theme(
#axis.title.x = element_blank(),
#axis.text.x = element_text(angle = 0, hjust=.5)
# legend.position = c(0.88, 0.83),
# legend.background=element_rect(fill = alpha("white", 0.1)),
legend.key.height = unit(0.7, "line") # Adjust this value to make the legend shorter
) +
scale_x_continuous(breaks = generate_vector(firstdate,lastdate),
minor_breaks = generate_vectorday(firstdate,lastdate),
limits = c(firstdate, lastdate) ) +
annotate("text",x = 51, y = weeky, label = "W08 ", size = 2.5) + #Week 08
annotate("text",x = 58, y = weeky, label = "W09 ", size = 2.5) + #Week 09
annotate("text",x = 65, y = weeky, label = "W10 ", size = 2.5) + #Week 10
annotate("text",x = 72, y = weeky, label = "W11 ", size = 2.5) + #Week 11
annotate("text",x = 79, y = weeky, label = "W12 ", size = 2.5) + #Week 12
annotate("text",x = 86, y = weeky, label = "W13 ", size = 2.5) + #Week 13
annotate("text",x = 93, y = weeky, label = "W14 ", size = 2.5) + #Week 14
annotate("text",x = 100, y = weeky, label = "W15 ", size = 2.5) + #Week 15
annotate("text",x = 107, y = weeky, label = "W16 ", size = 2.5) + #Week 16
annotate("text",x = 114, y = weeky, label = "W17 ", size = 2.5) + #Week 17
annotate("text",x = 121, y = weeky, label = "W18 ", size = 2.5) + #Week 18
annotate("text",x = 54, y = labely, label = "Steady state", size = 2.5) +
annotate("segment", y = segmenty, yend = segmenty, x = 47, xend = 62,
arrow = arrow(ends = "last", angle = 90, length = unit(.1,"cm"))) +
annotate("text",x = 65, y = labely, label = "Glycerol ", size = 2.5) +
annotate("segment", y = segmenty, yend = segmenty, x = 62, xend = 67,
arrow = arrow(ends = "both", angle = 90, length = unit(.1,"cm"))) +
annotate("text",x = 74, y = labely, label = " Inhibition", size = 2.5) +
annotate("segment", y = segmenty, yend = segmenty, x = 67, xend = 82,
arrow = arrow(ends = "both", angle = 90, length = unit(.1,"cm"))) +
annotate("text",x = 85, y = labely, label = " Recovery", size = 2.5) +
annotate("segment", y = segmenty, yend = segmenty, x = 82, xend = 90,
arrow = arrow(ends = "both", angle = 90, length = unit(.1,"cm"))) +
annotate("text",x = 94, y = labely, label = "Foaming ", size = 2.5) +
annotate("segment", y = segmenty, yend = segmenty, x = 90, xend = 98,
arrow = arrow(ends = "both", angle = 90, length = unit(.1,"cm"))) +
annotate("text",x = 101, y = labely, label = " Recovery", size = 2.5) +
annotate("segment", y = segmenty, yend = segmenty, x = 98, xend = 106,
arrow = arrow(ends = "both", angle = 90, length = unit(.1,"cm"))) +
annotate("text",x = 115, y = labely, label = "Steady state post foam", size = 2.5) +
annotate("segment", y = segmenty, yend = segmenty, x = 98, xend = 124,
arrow = arrow(ends = "first", angle = 90, length = unit(.1,"cm"))) +
annotate("text",x = 85, y = 4.25, label = "Feeding paused", size = 3, angle = 90, alpha = 0.5) +
annotate("rect", ymin = segmenty, ymax = 5, xmin = 82, xmax = 89, fill = "blue", alpha = .1) +
#annotate("text",x = as.Date("2023-06-03"), y = 301, label = "Foam", size = 2) +
# annotate("rect", ymin = 5.5, ymax = 125, xmin = as.Date("2023-06-01"), xmax = as.Date("2023-06-06"),
# fill = "blue", alpha = .1) +
ylab(expression(paste("Shannon"))) +
xlab("Day of operation")
ggsave("./Figures/Sup_lineplot_shannon.png", height=3, width=7.5, units='in', dpi=300)
Supplementary Figure S8
VFAs
weeky <- c(-150)
labely <- c(50)
segmenty <- c(0)
generate_vector <- function(x, y) {
result <- c()
while (x <= y) {
result <- c(result, x)
x <- x + 7
}
return(result)
}
generate_vectorday <- function(x, y) {
result <- c()
while (x <= y) {
result <- c(result, x)
x <- x + 1
}
return(result)
}
VFAs <- masterdata %>%
dplyr::select(SampleID.VFA, daysop, Date, AD, Treatment, Ethanol, Propanol,Butanol, X1.Hexanol, Acetic.acid, Propionic.acid, iso.Butyric.acid, Butyric.acid, iso.Valeric.acid, Valeric.acid, X4.Methyl.valeric.acid, Hexanoic.acid) %>% filter(SampleID.VFA != "")
# Control
p1 <- VFAs %>%
pivot_longer(Acetic.acid:Hexanoic.acid, names_to = "Compound") %>%
dplyr::filter(Treatment %in% c("Control")) %>%
# dplyr::filter(AD != c("AD3")) %>%
group_by(daysop, Compound, Treatment) %>%
summarise(value = median(value)) %>%
ggplot(., aes(x=daysop, y=value, fill=Compound)) +
geom_area(alpha=0.6 , linewidth=.1, colour="white", position = "stack") +
facet_wrap(vars(Treatment)) +
ylab(expression(paste("Concentration (mg L"^-1*")"))) +
scale_fill_viridis_d(labels = c("Acetic acid"," Butyric acid","Hexanoic acid","iso-Butyric acid","iso-Valeric acic","Propionic acid", "Valeric acid", "4-Methylvaleric acid")) +
theme_light() +
theme(
axis.title.x = element_blank(),
axis.text.x = element_text(angle = 0, hjust=1, vjust = 0.5)
) +
scale_x_continuous(breaks = generate_vector(firstdate,lastdate),
minor_breaks = generate_vectorday(firstdate,lastdate),
limits = c(firstdate, lastdate) ) +
annotate("text",x = 51, y = weeky, label = "W08 ", size = 2.5) + #Week 08
annotate("text",x = 58, y = weeky, label = "W09 ", size = 2.5) + #Week 09
annotate("text",x = 65, y = weeky, label = "W10 ", size = 2.5) + #Week 10
annotate("text",x = 72, y = weeky, label = "W11 ", size = 2.5) + #Week 11
annotate("text",x = 79, y = weeky, label = "W12 ", size = 2.5) + #Week 12
annotate("text",x = 86, y = weeky, label = "W13 ", size = 2.5) + #Week 13
annotate("text",x = 93, y = weeky, label = "W14 ", size = 2.5) + #Week 14
annotate("text",x = 100, y = weeky, label = "W15 ", size = 2.5) + #Week 15
annotate("text",x = 107, y = weeky, label = "W16 ", size = 2.5) + #Week 16
annotate("text",x = 114, y = weeky, label = "W17 ", size = 2.5) + #Week 17
annotate("text",x = 121, y = weeky, label = "W18 ", size = 2.5) + #Week 18
ylim(weeky, 5000) +
guides(fill = guide_legend(override.aes = list(size = 4),
keywidth = 0.5,
keyheight = 0.5)) +
annotate("text",x = 85, y = 3000, label = "Feeding paused", size = 3, angle = 90, alpha = .6) +
annotate("rect", ymin = 0, ymax = 5000, xmin = 82, xmax = 89, fill = "blue", alpha = .1)
# Treatment
p2 <- VFAs %>% pivot_longer(Acetic.acid:Hexanoic.acid, names_to = "Compound") %>%
dplyr::filter(Treatment %in% c("Treatment")) %>%
group_by(daysop, Compound, Treatment) %>%
summarise(value = median(value)) %>%
ggplot(., aes(x=daysop, y=value, fill=Compound)) +
geom_area(alpha=0.6 , linewidth=.1, colour="white", position = "stack") +
facet_wrap(vars(Treatment)) +
#ylab("Concentration (mg/L)") +
scale_fill_viridis_d(labels = c("Acetic acid"," Butyric acid","Hexanoic acid","iso-Butyric acid","iso-Valeric acic","Propionic acid", "Valeric acid", "4-Methylvaleric acid")) +
theme_light() +
theme(
axis.ticks.y = element_blank(),
axis.text.y = element_blank(),
axis.title.x = element_blank(),
axis.title.y = element_blank(),
axis.text.x = element_text(angle = 0, hjust=1, vjust = 0.5)) +
scale_x_continuous(breaks = generate_vector(firstdate,lastdate),
minor_breaks = generate_vectorday(firstdate,lastdate),
limits = c(firstdate, lastdate) ) +
annotate("text",x = 51, y = weeky, label = "W08 ", size = 2.5) + #Week 08
annotate("text",x = 58, y = weeky, label = "W09 ", size = 2.5) + #Week 09
annotate("text",x = 65, y = weeky, label = "W10 ", size = 2.5) + #Week 10
annotate("text",x = 72, y = weeky, label = "W11 ", size = 2.5) + #Week 11
annotate("text",x = 79, y = weeky, label = "W12 ", size = 2.5) + #Week 12
annotate("text",x = 86, y = weeky, label = "W13 ", size = 2.5) + #Week 13
annotate("text",x = 93, y = weeky, label = "W14 ", size = 2.5) + #Week 14
annotate("text",x = 100, y = weeky, label = "W15 ", size = 2.5) + #Week 15
annotate("text",x = 107, y = weeky, label = "W16 ", size = 2.5) + #Week 16
annotate("text",x = 114, y = weeky, label = "W17 ", size = 2.5) + #Week 17
annotate("text",x = 121, y = weeky, label = "W18 ", size = 2.5) + #Week 18
annotate("text",x = 64, y = 3000, label = "Glycerol added", size = 3, angle = 90, alpha = .6) +
annotate("rect",ymin =0, ymax = 5000, xmin = 62, xmax = 66, fill = "blue", alpha = .1) +
annotate("text",x = 85, y = 3000, label = "Feeding paused", size = 3, angle = 90, alpha = .6) +
annotate("rect", ymin = 0, ymax = 5000, xmin = 82, xmax = 89, fill = "blue", alpha = .1) +
annotate("text",x = 94, y = 3000, label = "Foam observed", size = 3, angle = 90, alpha = .6) +
annotate("rect", ymin = 0, ymax = 5000, xmin =92, xmax = 97, fill = "blue", alpha = .1) +
guides(fill = guide_legend(override.aes = list(size = 4),
keywidth = 0.5,
keyheight = 0.5))
g1 <- ggarrange(p1, theme_void(), p2, common.legend = TRUE, nrow = 1, legend = "top", align = "v",
widths = c(1,-0.02,0.9) )
g1 <- annotate_figure(g1, bottom = textGrob("Day of operation", gp = gpar(cex = 0.9)))
ggsave(plot = g1, "./Figures/areaplot_vfas_daysop.pdf", height=9, width=17, units='cm', dpi=300)
Supp. figure S9 CN
### CN
masterdata %>%
dplyr::filter(AD != "PSTWAS" &
AD != "Full-scale" &
week %in% c(9,10,11,12,13,14,15,16,17,18,19,20,21,22,23,24)) %>%
# drop_na(CN) %>%
#dplyr::filter(CN >= 0) %>%
ggboxplot(x = "week",
y = "CN",
fill = "Treatment",
legend = "right") +
ylab("C/N ratio") +
theme_light() +
scale_fill_manual("Treatment", values = cols) +
labs(x = "Week") +
geom_point(aes(x =week,y=CN, fill=Treatment),
position = position_jitterdodge(), alpha = 0.3) +
scale_x_discrete(breaks = c(9,10,11,12,13,14,15,16,17,18,21,24),
labels = c("W09","W10","W11","W12","W13","W14","W15","W16","W17","W18","W21","W24")) +
annotate("text",x = 5, y = 6, label = "Feeding\npaused", size = 2.5) +
theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))
ggsave( "./Figures/Sup_boxplot_CN.png", height=7, width=17,units='cm', dpi=600)
str(masterdata)
## 'data.frame': 1023 obs. of 68 variables:
## $ X : int 1 2 3 4 5 11 6 7 8 9 ...
## $ Date : Date, format: "2023-03-01" "2023-03-02" ...
## $ AD : Factor w/ 8 levels "AD1","AD2","AD3",..: 8 1 2 5 6 8 1 2 3 4 ...
## $ SludgeType : Factor w/ 5 levels "Control","Treatment",..: 5 1 2 2 2 5 1 2 1 1 ...
## $ Treatment : Factor w/ 4 levels "Control","Treatment",..: 4 1 2 2 2 4 1 2 1 1 ...
## $ daysop : int 0 1 1 1 1 2 2 2 2 2 ...
## $ week : int 1 1 1 1 1 1 1 1 1 1 ...
## $ Period : Factor w/ 8 levels "Converging","SteadyState",..: NA 1 1 1 1 NA 1 1 1 1 ...
## $ UQ.ID : chr "" "" "" "" ...
## $ Mtox.ID : chr "" "" "" "" ...
## $ EC50 : num NA NA NA NA NA NA NA NA NA NA ...
## $ SampleID.TSVS : chr "F7" "S13_AD1" "S13_AD1" "S14_AD5" ...
## $ TS_ww : num 3.35 1.09 1.76 1.82 2.36 ...
## $ VS_wwTS : num 86.1 66.7 68.3 67.6 67.8 ...
## $ VS_ww : num 2.89 0.73 1.21 1.24 1.6 ...
## $ VSr_pct : num NA NA NA NA NA NA NA NA NA NA ...
## $ SampleID.pH : chr "" "" "" "" ...
## $ pH : num NA NA NA NA NA NA NA NA NA NA ...
## $ SampleID.CODNH4 : chr "" "" "" "" ...
## $ COD_mgL : num NA NA NA NA NA NA NA NA NA NA ...
## $ TNH3_mgL : num NA NA NA NA NA NA NA NA NA NA ...
## $ FAN_mgL : num NA NA NA NA NA NA NA NA NA NA ...
## $ CODr_pct : num NA NA NA NA NA NA NA NA NA NA ...
## $ CH4_pct : int NA NA NA NA NA NA NA NA NA NA ...
## $ CH4_L : num NA NA NA NA NA NA NA NA NA NA ...
## $ CH4_L_kgVS : num NA NA NA NA NA NA NA NA NA NA ...
## $ O2_pct : num NA NA NA NA NA NA NA NA NA NA ...
## $ Gasflow_mL : num NA NA NA NA NA ...
## $ SampleID.CHN : chr "" "" "" "" ...
## $ C : num NA NA NA NA NA NA NA NA NA NA ...
## $ H_pct : num NA NA NA NA NA NA NA NA NA NA ...
## $ N_pct : num NA NA NA NA NA NA NA NA NA NA ...
## $ CN : num NA NA NA NA NA NA NA NA NA NA ...
## $ HC : num NA NA NA NA NA NA NA NA NA NA ...
## $ gC_kg : num NA NA NA NA NA NA NA NA NA NA ...
## $ gH_kg : num NA NA NA NA NA NA NA NA NA NA ...
## $ gN_kg : num NA NA NA NA NA NA NA NA NA NA ...
## $ SampleID.VFA : chr "" "" "" "" ...
## $ Ethanol : num NA NA NA NA NA NA NA NA NA NA ...
## $ Propanol : num NA NA NA NA NA NA NA NA NA NA ...
## $ Butanol : num NA NA NA NA NA NA NA NA NA NA ...
## $ X1.Hexanol : num NA NA NA NA NA NA NA NA NA NA ...
## $ Acetic.acid : num NA NA NA NA NA NA NA NA NA NA ...
## $ Propionic.acid : num NA NA NA NA NA NA NA NA NA NA ...
## $ iso.Butyric.acid : num NA NA NA NA NA NA NA NA NA NA ...
## $ Butyric.acid : num NA NA NA NA NA NA NA NA NA NA ...
## $ iso.Valeric.acid : num NA NA NA NA NA NA NA NA NA NA ...
## $ Valeric.acid : num NA NA NA NA NA NA NA NA NA NA ...
## $ X4.Methyl.valeric.acid: num NA NA NA NA NA NA NA NA NA NA ...
## $ Hexanoic.acid : num NA NA NA NA NA NA NA NA NA NA ...
## $ DNATubeID : chr "" "" "" "" ...
## $ SampleID.DNA : chr "" "" "" "" ...
## $ ngmL_qb : num NA NA NA NA NA ...
## $ DNA_mg_gCOD : num NA NA NA NA NA NA NA NA NA NA ...
## $ DNA_mg_gVS : num NA NA NA NA NA ...
## $ DNA_mg_L : chr "" "" "" "" ...
## $ DNA_mgLSludge : num NA NA NA NA NA 26.3 NA NA NA NA ...
## $ Observed : num NA NA NA NA NA 1130 NA NA NA NA ...
## $ Chao1 : num NA NA NA NA NA ...
## $ se.chao1 : num NA NA NA NA NA ...
## $ ACE : num NA NA NA NA NA ...
## $ se.ACE : num NA NA NA NA NA 16.8 NA NA NA NA ...
## $ Shannon : num NA NA NA NA NA 6.12 NA NA NA NA ...
## $ Simpson : num NA NA NA NA NA 1 NA NA NA NA ...
## $ InvSimpson : num NA NA NA NA NA ...
## $ Fisher : num NA NA NA NA NA ...
## $ pielou_ev : num NA NA NA NA NA 0.87 NA NA NA NA ...
## $ VFAS : num NA NA NA NA NA NA NA NA NA NA ...
Supp. figure S10 foam heights
firstdate <- 47 # selected periods for general manuscript
lastdate <- 124
weeky <- c(-8)
labely <- c(8)
segmenty <- c(0)
generate_vector <- function(x, y) {
result <- c()
while (x <= y) {
result <- c(result, x)
x <- x + 7
}
return(result)
}
generate_vectorday <- function(x, y) {
result <- c()
while (x <= y) {
result <- c(result, x)
x <- x + 1
}
return(result)
}
foamdata %>%
dplyr::filter(SludgeType %notin% c("Full-scale")) %>% # filter columns out
ggline(x = "DayofOp",
y = "Max",
color = "SludgeType",
legend = "right",
add = "mean_se",
add.params = list(width = 1)) +
# scale_color_manual(values = cols) + # create a cols vector for custom colours
theme_light() +
scale_color_manual("Sludge Supernatent", values = cols) +
theme(
axis.text.x = element_text(angle = 0, hjust=0.5),
legend.position = c(0.90, 0.825),
legend.background=element_rect(fill = alpha("white", 0.5)),
legend.ticks.length = element_blank(),
legend.title = element_blank() ) +
scale_x_continuous(breaks = generate_vector(firstdate,lastdate),
minor_breaks = generate_vectorday(firstdate,lastdate),
limits = c(firstdate, lastdate) ) +
annotate("text",x = 51, y = weeky, label = "W08 ", size = 2.5) + #Week 08
annotate("text",x = 58, y = weeky, label = "W09 ", size = 2.5) + #Week 09
annotate("text",x = 65, y = weeky, label = "W10 ", size = 2.5) + #Week 10
annotate("text",x = 72, y = weeky, label = "W11 ", size = 2.5) + #Week 11
annotate("text",x = 79, y = weeky, label = "W12 ", size = 2.5) + #Week 12
annotate("text",x = 86, y = weeky, label = "W13 ", size = 2.5) + #Week 13
annotate("text",x = 93, y = weeky, label = "W14 ", size = 2.5) + #Week 14
annotate("text",x = 100, y = weeky, label = "W15 ", size = 2.5) + #Week 15
annotate("text",x = 107, y = weeky, label = "W16 ", size = 2.5) + #Week 16
annotate("text",x = 114, y = weeky, label = "W17 ", size = 2.5) + #Week 17
annotate("text",x = 121, y = weeky, label = "W18 ", size = 2.5) + #Week 18
annotate("text",x = 54, y = labely, label = "Steady state", size = 2.5) +
annotate("segment", y = segmenty, yend = segmenty, x = 47, xend = 62,
arrow = arrow(ends = "last", angle = 90, length = unit(.1,"cm"))) +
# annotate("text",x = 65, y = labely, label = "Glycerol ", size = 2.5) +
annotate("segment", y = segmenty, yend = segmenty, x = 62, xend = 67,
arrow = arrow(ends = "both", angle = 90, length = unit(.1,"cm"))) +
annotate("text",x = 74, y = labely, label = " Inhibition", size = 2.5) +
annotate("segment", y = segmenty, yend = segmenty, x = 67, xend = 82,
arrow = arrow(ends = "both", angle = 90, length = unit(.1,"cm"))) +
annotate("text",x = 86, y = labely, label = "Recovery", size = 2.5) +
annotate("segment", y = segmenty, yend = segmenty, x = 82, xend = 90,
arrow = arrow(ends = "both", angle = 90, length = unit(.1,"cm"))) +
#annotate("text",x = 94, y = labely, label = "Foaming ", size = 2.5) +
annotate("segment", y = segmenty, yend = segmenty, x = 90, xend = 98,
arrow = arrow(ends = "both", angle = 90, length = unit(.1,"cm"))) +
annotate("text",x = 101, y = labely, label = " Recovery", size = 2.5) +
annotate("segment", y = segmenty, yend = segmenty, x = 98, xend = 106,
arrow = arrow(ends = "both", angle = 90, length = unit(.1,"cm"))) +
annotate("text",x = 115, y = labely, label = "Steady state post foam", size = 2.5) +
annotate("segment", y = segmenty, yend = segmenty, x = 98, xend = 124,
arrow = arrow(ends = "first", angle = 90, length = unit(.1,"cm"))) +
annotate("text",x = 85, y = 100, label = "Feeding paused", size = 4, angle = 90, alpha = 0.5) +
annotate("rect", ymin = 0, ymax = 200, xmin = 82, xmax = 89, fill = "blue", alpha = .1) +
annotate("text",x = 64, y = 75, label = "Glycerol", size = 4, angle = 90, , alpha = 0.5) +
annotate("rect",ymin =0.5, ymax = 150, xmin = 62, xmax = 67,
fill = "blue", alpha = .1) +
annotate("text",x = 93.5, y = 75, label = "Foam observed", size = 3.5, angle = 90, , alpha = 0.5) +
annotate("rect", ymin = 0.5, ymax = 150, xmin = 91, xmax = 96,
fill = "blue", alpha = .1) +
ylab("Max foam height (mm)") +
xlab("Day of operation")
ggsave("./Figures/Sup_lineplot_foam_daysop.png", height=2.5, width=6, units='in', dpi=300)
Supp. figure S12 phylum abund
Control
firstdate <- -2 # selected periods for Supplementary
lastdate <- 166
weeky <- c(-5)
labely <- c(.05)
segmenty <- c(0)
generate_vector <- function(x, y) {
result <- c()
while (x <= y) {
result <- c(result, x)
x <- x + 7
}
return(result)
}
generate_vectorday <- function(x, y) {
result <- c()
while (x <= y) {
result <- c(result, x)
x <- x + 1
}
return(result)
}
ps.rel <- microbiome::aggregate_taxa(ps.flt, "Phylum")
ps.rel <- ps.rel %>%
microbiome::transform(transform = 'compositional')
df <- phyloseq::psmelt(ps.rel) %>%
mutate(Abundance = Abundance *100)
p1 <- df %>%
group_by(daysop, Kingdom, Phylum, CSTR, SludgeType) %>%
summarise(value = median(Abundance)) %>%
dplyr::filter(SludgeType %in% c("Control")) %>%
ggplot(., aes(x=daysop, y=value, fill=Phylum)) +
geom_area(alpha=0.6 , linewidth=.1, colour="white", position = "stack") +
facet_wrap(vars(CSTR), ncol = 1) +
scale_x_continuous(breaks = generate_vector(firstdate,lastdate),
minor_breaks = generate_vectorday(firstdate,lastdate),
limits = c(firstdate, lastdate) ) +
ylab("Relative Abundance") +
theme_light() +
# scale_color_manual("Treatment", values = cols) +
theme(
axis.title.x = element_blank(),
#axis.text.x = element_blank(),
axis.text.x = element_text(angle = 90, hjust=1),
legend.position = "top",
legend.title = element_blank()
) +
# annotate("text",x = 2, y = weeky, label = "W01 ", size = 2.5) + #Week 01
annotate("text",x = 9, y = weeky, label = "W02 ", size = 2.5) + #Week 02
# annotate("text",x = 16, y = weeky, label = "W03 ", size = 2.5) + #Week 03
annotate("text",x = 23, y = weeky, label = "W04 ", size = 2.5) + #Week 04
# annotate("text",x = 30, y = weeky, label = "W05 ", size = 2.5) + #Week 05
annotate("text",x = 37, y = weeky, label = "W06 ", size = 2.5) + #Week 06
# annotate("text",x = 44, y = weeky, label = "W07 ", size = 2.5) + #Week 07
annotate("text",x = 51, y = weeky, label = "W08 ", size = 2.5) + #Week 08
# annotate("text",x = 58, y = weeky, label = "W09 ", size = 2.5) + #Week 09
annotate("text",x = 65, y = weeky, label = "W10 ", size = 2.5) + #Week 10
# annotate("text",x = 72, y = weeky, label = "W11 ", size = 2.5) + #Week 11
annotate("text",x = 79, y = weeky, label = "W12 ", size = 2.5) + #Week 12
# annotate("text",x = 86, y = weeky, label = "W13 ", size = 2.5) + #Week 13
annotate("text",x = 93, y = weeky, label = "W14 ", size = 2.5) + #Week 14
# annotate("text",x = 100, y = weeky, label = "W15 ", size = 2.5) + #Week 15
annotate("text",x = 107, y = weeky, label = "W16 ", size = 2.5) + #Week 16
# annotate("text",x = 114, y = weeky, label = "W17 ", size = 2.5) + #Week 17
annotate("text",x = 121, y = weeky, label = "W18 ", size = 2.5) + #Week 18
# annotate("text",x = 128, y = weeky, label = "W19 ", size = 2.5) + #Week 19
annotate("text",x = 135, y = weeky, label = "W20 ", size = 2.5) + #Week 20
# annotate("text",x = 142, y = weeky, label = "W21 ", size = 2.5) + #Week 21
annotate("text",x = 149, y = weeky, label = "W22 ", size = 2.5) + #Week 22
# annotate("text",x = 156, y = weeky, label = "W23 ", size = 2.5) + #Week 23
annotate("text",x = 163, y = weeky, label = "W24 ", size = 2.5) + #Week 24
annotate("text",x = 85, y = 50, label = "Feeding paused", size = 4, angle = 90, alpha = 0.6) +
annotate("rect", ymin = 0, ymax = 100, xmin = 82, xmax = 89,
fill = "blue", alpha = .1) +
guides(fill=guide_legend(ncol=6)) +
ggtitle("Control")
Treatment
ps.rel <- microbiome::aggregate_taxa(ps.flt, "Phylum")
ps.rel <- ps.rel %>%
microbiome::transform(transform = 'compositional')
df <- phyloseq::psmelt(ps.rel) %>%
mutate(Abundance = Abundance *100)
p2 <- df %>%
group_by(daysop, Kingdom, Phylum, CSTR, SludgeType) %>%
summarise(value = median(Abundance)) %>%
dplyr::filter(SludgeType %in% c("Treatment")) %>%
ggplot(., aes(x=daysop, y=value, fill=Phylum)) +
geom_area(alpha=0.6 , linewidth=.1, colour="white", position = "stack") +
facet_wrap(vars(CSTR), ncol = 1) +
scale_x_continuous(breaks = generate_vector(firstdate,lastdate),
minor_breaks = generate_vectorday(firstdate,lastdate),
limits = c(firstdate, lastdate) ) +
ylab("Relative Abundance") +
theme_light() +
theme(
axis.title.x = element_blank(),
axis.title.y = element_blank(),
axis.text.y = element_blank(),
axis.text.x = element_text(angle = 90, hjust=1),
legend.position = "top",
legend.title = element_blank()
) +
# annotate("text",x = 2, y = weeky, label = "W01 ", size = 2.5) + #Week 01
annotate("text",x = 9, y = weeky, label = "W02 ", size = 2.5) + #Week 02
# annotate("text",x = 16, y = weeky, label = "W03 ", size = 2.5) + #Week 03
annotate("text",x = 23, y = weeky, label = "W04 ", size = 2.5) + #Week 04
# annotate("text",x = 30, y = weeky, label = "W05 ", size = 2.5) + #Week 05
annotate("text",x = 37, y = weeky, label = "W06 ", size = 2.5) + #Week 06
# annotate("text",x = 44, y = weeky, label = "W07 ", size = 2.5) + #Week 07
annotate("text",x = 51, y = weeky, label = "W08 ", size = 2.5) + #Week 08
# annotate("text",x = 58, y = weeky, label = "W09 ", size = 2.5) + #Week 09
annotate("text",x = 65, y = weeky, label = "W10 ", size = 2.5) + #Week 10
# annotate("text",x = 72, y = weeky, label = "W11 ", size = 2.5) + #Week 11
annotate("text",x = 79, y = weeky, label = "W12 ", size = 2.5) + #Week 12
# annotate("text",x = 86, y = weeky, label = "W13 ", size = 2.5) + #Week 13
annotate("text",x = 93, y = weeky, label = "W14 ", size = 2.5) + #Week 14
# annotate("text",x = 100, y = weeky, label = "W15 ", size = 2.5) + #Week 15
annotate("text",x = 107, y = weeky, label = "W16 ", size = 2.5) + #Week 16
# annotate("text",x = 114, y = weeky, label = "W17 ", size = 2.5) + #Week 17
annotate("text",x = 121, y = weeky, label = "W18 ", size = 2.5) + #Week 18
# annotate("text",x = 128, y = weeky, label = "W19 ", size = 2.5) + #Week 19
annotate("text",x = 135, y = weeky, label = "W20 ", size = 2.5) + #Week 20
# annotate("text",x = 142, y = weeky, label = "W21 ", size = 2.5) + #Week 21
annotate("text",x = 149, y = weeky, label = "W22 ", size = 2.5) + #Week 22
# annotate("text",x = 156, y = weeky, label = "W23 ", size = 2.5) + #Week 23
annotate("text",x = 163, y = weeky, label = "W24 ", size = 2.5) + #Week 24
annotate("text",x = 85, y = 50, label = "Feeding paused", size = 4, angle = 90, alpha = 0.6) +
annotate("rect", ymin = 0, ymax = 100, xmin = 82, xmax = 89,
fill = "blue", alpha = .1) +
annotate("text",x = 64, y = 50, label = "Glycerol", size = 4, angle = 90, alpha = 0.6) +
annotate("rect",ymin =0, ymax = 100, xmin = 62, xmax = 67,
fill = "blue", alpha = .1) +
annotate("text",x = 94, y = 50, label = "Foam observed", size = 4, angle = 90, alpha = 0.65) +
annotate("rect",ymin =0, ymax = 100, xmin = 93, xmax = 96,
fill = "blue", alpha = .1) +
guides(fill=guide_legend(ncol=3)) +
ggtitle("Treatment")
Combined
g1 <- ggarrange(p1, p2, common.legend = TRUE, widths = c(0.52,0.48), legend = "top")
g1 <- annotate_figure(g1, bottom = textGrob("Day of operation", gp = gpar(cex = 0.9)))
ggsave( "./Figures/Sup_areaplot_all_phylum.pdf", height=27, width=26, units='cm', dpi=300)
Supp. figure S13 Ampviz heatmaps 16S top30
psobject <- ps.flt
set.seed(1000)
psobject <- phyloseq::rarefy_even_depth(psobject)
psobject <- phyloseq::prune_samples(sample_data(psobject)$Date %in% c("2023-03-03","2023-03-15", "2023-04-12", "2023-04-26","2023-05-01", "2023-05-15", "2023-05-24", "2023-06-02","2023-06-09","2023-06-05","2023-06-14","2023-06-28","2023-07-26","2023-08-09"), psobject)
psobject <- phyloseq::prune_samples(sample_data(psobject)$SampleID.TSVS %notin% c("F16", "F77", "F205", "F150"), psobject)
psobject <- phyloseq::prune_samples(sample_data(psobject)$DNATubeID %notin% c("#061", "#062","#063", "#064","#065","#066",
"#216", "#217","#218", "#219","#220","#221"), psobject)
metadata <- data.frame(sample_data(psobject)) %>%
rownames_to_column("#OTU ID")
taxtable <- data.frame(tax_table(psobject))
ASVtable <- data.frame(otu_table(psobject))
amp <- amp_load(otutable = ASVtable,
metadata = metadata,
taxonomy = taxtable)
p2 <- amp_heatmap(amp, tax_aggregate = "Phylum", tax_show = 30, showRemainingTaxa = TRUE,
normalise = TRUE, facet_by = "SludgeType", group_by = "daysop")
p2 <- p2 + theme(
strip.text.x = element_text(size = 8, colour = "black", angle = 90),
axis.text.x = element_text(size = 12)
)
ggsave(plot = p2 , "./Figures/Sup_heatmap_v3v4midas53_all_phylum.png", height=18, width=29, units='cm', dpi=200)
Supp. figure S14 Ampviz heatmaps metagenome phyla
library(ampvis2)
# Squeezemeta as described
psobject <- psmg
sample_data(psobject)$SludgeType <- factor(sample_data(psobject)$SludgeType, levels = c("Control", "Treatment", "Foam", "Full-scale", "PSTWAS"))
#set.seed(1000)
#psobject <- phyloseq::rarefy_even_depth(psobject)
metadata <- data.frame(sample_data(psobject)) %>%
rownames_to_column("#OTU ID")
#metadata$`#OTU ID` <- str_replace_all(metadata$`#OTU ID`, "-", ".")
taxtable <- data.frame(tax_table(psobject)) %>% dplyr::rename("Kingdom" = SuperKingdom)
ASVtable <- data.frame(otu_table(psobject))
amp <- amp_load(otutable = ASVtable,
metadata = metadata,
taxonomy = taxtable)
p <- amp_heatmap(amp, tax_aggregate = "Phylum", tax_show = 30, showRemainingTaxa = TRUE,
normalise = TRUE, facet_by = "SludgeType", group_by = "daysop")
pmg <- p + theme(strip.text.x = element_text(size = 8, colour = "black", angle = 90))
## HMMER Predicted, 16S extracted, Midas 53
df <- read.csv("./data/2024-06-11_arcbac_phylotabel.csv")
df <- df %>%
dplyr::rename(S561toS566_AD.AD6 = Chris.11,
S1001_3_4_AD1_3_4 = Chris.12,
S1002_5_6_AD2_5_6 = Chris.13,
S1147_49_50_AD1_3_4 = Chris.14,
S1148_51_52_AD2_5_6 = Chris.15,
S1154_57_58_AD2_5_6 = Chris.16,
F5 = Chris.01,
F55 = Chris.02,
F110 = Chris.03,
F124 = Chris.04,
F179 = Chris.05,
S119.full = Chris.06,
S370.full = Chris.07,
S1273.ADFull = Chris.08,
S1554.ADFull = Chris.09,
S1813.ADFull = Chris.10)
ncol(df %>% select( F5:S1154_57_58_AD2_5_6) )
## [1] 16
colSums(df %>% select( F5:S1154_57_58_AD2_5_6) )
## F5 F55 F110 F124
## 12436 11891 11201 11643
## F179 S119.full S370.full S1273.ADFull
## 11173 14465 12257 12843
## S1554.ADFull S1813.ADFull S561toS566_AD.AD6 S1001_3_4_AD1_3_4
## 12203 12506 13281 11726
## S1002_5_6_AD2_5_6 S1147_49_50_AD1_3_4 S1148_51_52_AD2_5_6 S1154_57_58_AD2_5_6
## 14219 13718 14385 14774
amp <- amp_load(otutable = df,
metadata = (metadata) )
p <- amp_heatmap(amp, tax_aggregate = "Phylum", tax_show = 30, showRemainingTaxa = TRUE,
normalise = TRUE, group_by = "daysop", facet_by = "SludgeType") #+ ggtitle("HMMER/sintax classification")
phmmer.midas <- p + theme(strip.text.x = element_text(size = 8, colour = "black", angle = 90))
g1 <- ggarrange(pmg, phmmer.midas, ncol = 1, align = "hv", labels = "AUTO")
ggsave(plot = g1, "./Figures/Sup_heatmap_hmmer_squemeta_phylum_comparison.png", height=33, width=25, units='cm', dpi=200)
Supp. figure S15 Picrust2 heatmap
library(circlize)
library(ComplexHeatmap)
output <- readRDS("./data/ancombcoutputpicrust-09Jul24")
res_prim = output$res
filtervec <- (res_prim %>% dplyr::filter(q_TreatmentTreatment <= 0.05))$taxon
ps.tmp <- ps_pw
ps.tmp <- prune_taxa(taxa_sums(ps.tmp) > 10000, ps.tmp) # minimum reads per feature
ps.tmp <- subset_samples(ps.tmp, sample_data(ps.tmp)$Period %in% c("Inhibition","Recovery/Feedingpause", "Recovery/Foaming"))
ps.tmp <- subset_samples(ps.tmp, sample_data(ps.tmp)$SludgeType %in% c("Control", "Treatment"))
ps.tmp <- prune_taxa(taxa_sums(ps.tmp) > 0, ps.tmp) # minimum reads per feature
ps.tmp <- microbiome::transform(ps.tmp, transform = "hellinger")
ps.tmp <- phyloseq::subset_taxa(ps.tmp, taxa_names(ps.tmp) %in% filtervec)
meta.tmp <- sample_data(ps.tmp)
sample_data(ps.tmp)$CSTR <- str_replace_all(sample_data(ps.tmp)$AD, c("AD1"="C1",
"AD2"="T1",
"AD3"="C2",
"AD4"="C3",
"AD5"="T2",
"AD6"="T3",
"AD_full" = "Full-scale") )
vars <- c("CSTR", "Treatment")
sample_data(ps.tmp) <- tidyr::unite(data.frame(sample_data(ps.tmp)), "CSTR_Treat", all_of(vars),remove = FALSE )
summary(sample_data(ps.tmp)$Treatment)
## Control Treatment
## 27 27
summary(sample_data(ps.tmp)$AD)
## AD1 AD2 AD3 AD4 AD5 AD6
## 9 9 9 9 9 9
# plot
taxtable <- data.frame(tax_table(ps.tmp))
colAnn <- HeatmapAnnotation(
`Gasflow (ml/day)` = anno_barplot(sample_data(ps.tmp)$Gasflow_mL, width = unit(4, "cm")),
Digester = as.character(sample_data(ps.tmp)$Treatment),
col = list(Digester= c( "Control" = "#4D98AC", "Treatment" = "#985C64")
), annotation_name_gp= gpar(fontsize = 8), show_legend = FALSE
)
set.seed(1232145)
row_ha <- HeatmapAnnotation(
which = c("row"),
`MetaCyc class` = taxtable$CLASS1, annotation_name_gp= gpar(fontsize = 8),
annotation_legend_param = list(`MetaCyc class` = list(ncol = 1)),
show_legend = TRUE
)
fig <- Heatmap((otu_table(ps.tmp)),
name = "Abundances", #title of legend
row_names_gp = gpar(fontsize = 6.5),
#split = sample_data(ps.tmp)$Treatment,
#col_fun,
show_heatmap_legend = FALSE,
show_column_names = TRUE,
show_row_names = TRUE,
row_labels = c(data.frame(tax_table(ps.tmp))$Pathway2),
column_labels = c(data.frame(sample_data(ps.tmp))$daysop),
column_names_gp = gpar(fontsize = 6),
cluster_columns = TRUE,
clustering_method_columns = "ward.D",
top_annotation = colAnn,
left_annotation = row_ha
#row_km = 1, column_split = 2
)
#column_split = 4
# Text size for row names
png("./Figures/Sup_heatmap_metacyc_picrust_noleg.png", width = 19, height = 20,pointsize = 12,units="cm",res=300)
draw(fig, padding = unit(c(3, 3, 2, 7), "mm"), merge_legend = TRUE, heatmap_legend_side = c("top"),
annotation_legend_side = "bottom") ## see right heatmap in following
Supp. figure S16 - Metabolome PCA - is done in
EcologyofFoaming_figure7
Supp. figure S18 Kegg pathways
Aldex and create filtervec
library(ALDEx2)
# Filter out full-scale AD
taxtable <- data.frame(tax_table(ps_kegg))
sample_data(ps_kegg)$SludgeType3 <- str_replace_all(sample_data(ps_kegg)$SludgeType,
c("Control"="RMITAD",
"Treatment"="RMITAD",
"Foam"="RMITAD",
"PSTWAS"="PSTWAS",
"Full-scale"="Full-scale") )
sample_data(ps_kegg)$SludgeType4 <- str_replace_all(sample_data(ps_kegg)$SludgeType,
c("Control"="Control",
"Treatment"="Treatment",
"Foam"="Treatment",
"PSTWAS"="PSTWAS",
"Full-scale"="Full-scale") )
ps.tmp <- subset_samples(ps_kegg, sample_data(ps_kegg)$SludgeType3 %in% c("RMITAD"))
ps.tmp = subset_taxa(ps.tmp, PathwayL1 %notin% c("Human Diseases", "Organismal Systems"))
ps.tmp <- prune_taxa(taxa_sums(ps.tmp) >= 100, ps.tmp) # minimum reads per feature
prev.kegg <- prevalencedf(ps.tmp , Pathway)
abun <- data.frame(phyloseq::otu_table(ps.tmp))
conds <- as.character(sample_data(ps.tmp)$SludgeType4)
summary(as.factor(conds))
## Control Treatment
## 3 3
x <- aldex.clr(abun, conds, mc.samples=1000, denom="all",
verbose = TRUE)
x.tt <- aldex.ttest(x, paired.test=FALSE, verbose=FALSE)
x.effect <- aldex.effect(x, CI=T, verbose=FALSE)
x.all <- data.frame(x.tt,x.effect) %>% rownames_to_column("ko")
png("./Figures/Sup_aldex2.png", width = 14, height = 12,pointsize = 12,units="cm",res=300)
aldex.plot(x.all, type="MA", test="wilcox")
dev.off()
## quartz_off_screen
## 2
#sgn <- sign(x.effect$effect.low) == sign(x.effect$effect.high)
#par(mfrow=c(1,2))
#plot(x.effect$rab.all, x.effect$diff.btw, pch=19, cex=0.3, col="grey", xlab="Abundance", ylab="Difference", main="Bland-Altman")
#points(x.effect$rab.all[abs(x.effect$effect) >=2], x.effect$diff.btw[abs(x.effect$effect) >=2], pch=19, cex=0.5, col="red")
#points(x.effect$rab.all[sgn], x.effect$diff.btw[sgn], cex=0.7, col="blue")
#plot(x.effect$diff.win, x.effect$diff.btw, pch=19, cex=0.3, col="grey",xlab="Dispersion", ylab="Difference", main="Effect")
#points(x.effect$diff.win[abs(x.effect$effect) >=2], x.effect$diff.btw[abs(x.effect$effect) >=2], pch=19, cex=0.5, col="red")
#points(x.effect$diff.win[sgn], x.effect$diff.btw[sgn], cex=0.7, col="blue")
#abline(0,2, lty=2, col="grey")
#abline(0,-2, lty=2, col="grey")
# filtervec
# by effect size
# filtervec <- (x.all %>% dplyr::filter(effect >= 2 |
# effect <= -2))
# by p value
filtervec <- (x.all %>% dplyr::filter(wi.eBH <= .05))
filtervec <- filtervec$ko
Complex Heatmap
library(circlize)
library(ComplexHeatmap)
prev.kegg <- prevalencedf(ps_kegg, Pathway)
sample_data(ps_kegg)$SludgeType3 <- str_replace_all(sample_data(ps_kegg)$SludgeType,
c("Control"="RMITAD",
"Treatment"="RMITAD",
"Foam"="RMITAD",
"PSTWAS"="PSTWAS",
"Full-scale"="Full-scale") )
sample_data(ps_kegg)$SludgeType3<- factor(sample_data(ps_kegg)$SludgeType3)
# filter
ps.tmp <- subset_samples(ps_kegg, sample_data(ps_kegg)$SludgeType3 %in% c("RMITAD"))
ps.tmp <- microbiome::transform(ps.tmp, transform = "clr")
ps.tmp <- phyloseq::subset_taxa(ps.tmp, taxa_names(ps.tmp) %in% filtervec)
#ps.tmp <- prune_taxa(taxa_sums(ps.tmp) > 0, ps.tmp) # minimum reads per feature
# plot
#col_fun = colorRamp2(c(min(otu_table(ps.tmp)), mean(otu_table(ps.tmp)), max(otu_table(ps.tmp))), c("black", "steelblue", "#FCFDBFFF"))
taxtable <- data.frame(tax_table(ps.tmp))
colAnn <- HeatmapAnnotation(
#Date = as.character(sample_data(ps.tmp)$Date),
`Day of operation` = as.character(sample_data(ps.tmp)$daysop),
col = list(`Day of operation` = c( "61" = "#D2A277","75" = "#E495A5",
"93" = "#72BB83"))
)
set.seed(1232145)
row_ha = rowAnnotation(`KEGG pathway` = taxtable$PathwayL2,
annotation_name_gp= gpar(fontsize = 8))
#cols <- colorspace::rainbow_hcl(length(unique(sample_data(ps.tmp)$Date)))
#row_ano = rowAnnotation(foo = anno_block(gp = gpar(fill = 2:4)))
fig <- Heatmap((otu_table(ps.tmp)),
name = "Abundances", #title of legend
row_names_gp = gpar(fontsize = 6.5),
#split = sample_data(ps.tmp)$Treatment,
#col_fun,
show_column_names = TRUE,
show_row_names = TRUE,
row_labels = c(data.frame(tax_table(ps.tmp))$Pathway),
column_labels = c(data.frame(sample_data(ps.tmp))$SludgeType),
column_names_gp = gpar(fontsize = 6.5),
clustering_method_columns = "ward.D",
top_annotation = colAnn,
left_annotation = row_ha,
row_km = 1
)
## Warning in class(matrix) <- "matrix": Setting class(x) to "matrix" sets
## attribute to NULL; result will no longer be an S4 object
#column_split = 4
# Text size for row names
fig
png("./Figures/Sup_heatmap_kegg_metagenome.png", width = 17, height = 12,pointsize = 12,units="cm",res=300)
draw(fig, padding = unit(c(3, 3, 2, 2), "mm"), merge_legend = TRUE) ## see right heatmap in following
dev.off()
## quartz_off_screen
## 2
# 550 x 850
Supp. figure S19 violin/heatmap filaments 16S
Violin
ps.temp <- ps.flt
set.seed(1234)
# summary(sample_sums(ps.temp))
# sort(sample_sums(ps.temp))
ps.temp <- phyloseq::rarefy_even_depth(ps.temp, sample.size = 30563, replace = FALSE)
#ps.temp <- ps.temp %>% microbiome::transform("clr")
#ps.temp <- microbiome::aggregate_taxa(ps.temp, level = "Genus") # faster
ps.temp <- ps.temp %>% microbiome::transform("compositional")
ps.summary <- ps.temp %>% phyloseq::psmelt() %>%
dplyr::select(Date, SludgeType, week, Treatment, Abundance, Kingdom, Phylum, Class, Order, Family, Genus, Species) %>%
dplyr::filter(Abundance > 0) %>%
group_by(Date,SludgeType, week, Treatment, Kingdom, Phylum, Class, Order, Family,Genus, Species) %>%
summarise(Abundance = mean(Abundance) *100) %>% unique()
ps.summary_sel <- ps.summary %>%
dplyr::filter( str_detect(Family, "Nocard") | str_detect(Family, "Microtrichaceae") | str_detect(Family, "Corynebacter") |
str_detect(Family, "Mycobacter") | str_detect(Family, "Tsukamurellaceae") |
str_detect(Genus, "Tetrasphaera") ) %>%
dplyr::filter(week %in% c(10, 11,12, 13, 14,15,16))
p <- ggpubr::ggviolin(ps.summary_sel, x = "week",
y = "Abundance",
ylab = "Abundance (%)",
xlab = "Week",
# facet.by = c("PMA_treat"),
fill = "SludgeType") +
# shape = "PMA_treat") +
theme_bw() +
scale_fill_manual("Sample type",values = cols, labels = c("Control CSTR", "Treatment CSTR", "Foam", "Full-scale CSTR", "PS/TWAS")) +
theme( panel.grid.major = element_blank(), panel.grid.minor = element_blank() ) +
geom_point(aes(x =week,y=Abundance, fill=SludgeType),
position = position_jitterdodge(), alpha = 0.15) +
theme(axis.text.x = element_text(angle = 0, hjust = 0.5, size = 12)) +
ggpubr::stat_compare_means(aes(x =week,y=Abundance, group=SludgeType,
label = paste0(..p.signif..)),
method = "kruskal")
heatmap
### heatmap of these
ps.summary <- ps.temp %>% phyloseq::psmelt() %>%
dplyr::select(OTU, Date, SludgeType, week, Treatment, Abundance, Kingdom, Phylum, Class, Order, Family, Genus, Species) %>%
dplyr::filter(Abundance > 0)
#ps.summary_sel <- ps.summary %>%
# dplyr::filter( str_detect(Family, "Nocard") | str_detect(Family, "Microtrichaceae") | str_detect(Family, "Corynebacter") |
# str_detect(Family, "Mycobacter") | str_detect(Family, "Tsukamurellaceae") ) %>%
# dplyr::filter(Period %in% c("SteadyState", "Inhibition", "Recovery/Foaming", "Recovery/Feedingpause",
# "Glycerol"))
ps.summary_sel <- ps.summary %>%
dplyr::filter( str_detect(Family, "Nocard") | str_detect(Family, "Microtrichaceae") | str_detect(Family, "Corynebacter") |
str_detect(Family, "Mycobacter") | str_detect(Family, "Tsukamurellaceae") |
str_detect(Genus, "Tetrasphaera") | str_detect(Genus, "Methanospirillum"))
library(ampvis2)
ps.temp <- ps.flt
psobject <- ps.temp
psobject <- subset_samples(psobject, sample_data(psobject)$week %in% c(10, 11,12, 13,14,15,16))
psobject <- prune_taxa(taxa_sums(psobject) > 0, psobject) #minimum reads per feature
taxtable <- data.frame(tax_table(psobject))
ASVtable <- data.frame(otu_table(psobject))
metaobject <- data.frame(sample_data(psobject)) %>%
rownames_to_column("#OTU ID")
amp <- amp_load(otutable = ASVtable,
metadata = metaobject,
taxonomy = taxtable)
## Warning: Could not find a column named one of "OTU, ASV, #OTU ID" in otutable.
## Using row names as OTU ID's.
## Warning: Could not find a column named one of "OTU, ASV, #OTU ID" in taxonomy.
## Using row names as OTU ID's.
filtervec <- ps.summary_sel$Genus %>% unique()
amp <- amp %>% amp_subset_taxa( normalise = TRUE, tax_vector = filtervec )
## 4385 OTUs have been filtered
## Before: 4450 OTUs
## After: 65 OTUs
p1 <- amp_heatmap(amp, tax_aggregate = "Genus", tax_add = c("Family"),
group_by = "week",
tax_show = 30, facet_by = "SludgeType",
showRemainingTaxa = TRUE, normalise = FALSE,
min_abundance = 0.001) +
theme(axis.text.x = element_text(size = 12))
## Warning: There are only 19 taxa, showing all
g <- ggarrange(p, p1, nrow = 2, labels = "AUTO")
## Warning in scale_fill_gradientn(colours = color.pal, trans = plot_colorscale, :
## log-10 transformation introduced infinite values.
g
ggsave(plot = g, "./Figures/Sup_violin_heatmap_filamen_hydrop.png", height=7, width=12, units='in', dpi=300)
Supp. figure S20 heatmap filaments squeezemeta
### ampviz heatmap
ps.temp <- psmg
sample_sums(ps.temp)
## F5 F55 F110 F124
## 39179916 34202748 35175590 35847078
## F179 S119.full S370.full S1273.ADFull
## 35462072 43047926 36930810 41085132
## S1554.ADFull S1813.ADFull S561toS566_AD.AD6 S1001_3_4_AD1_3_4
## 38015368 36078606 33312012 30583940
## S1002_5_6_AD2_5_6 S1147_49_50_AD1_3_4 S1148_51_52_AD2_5_6 S1154_57_58_AD2_5_6
## 35009536 34570590 32645912 35602780
sample_data(ps.temp)$SludgeType <- factor(sample_data(ps.temp)$SludgeType , levels = c("Control","Treatment","Foam","Full-scale", "PSTWAS"))
ps.summary <- ps.temp %>% phyloseq::psmelt() %>%
dplyr::select(OTU, Date, SludgeType, week, Abundance, SuperKingdom, Phylum, Class, Order, Family, Genus, Species) %>%
dplyr::filter(Abundance > 0)
#ps.summary_sel <- ps.summary %>%
# dplyr::filter( str_detect(Family, "Nocard") | str_detect(Family, "Microtrichaceae") | str_detect(Family, "Corynebacter") |
# str_detect(Family, "Mycobacter") | str_detect(Family, "Tsukamurellaceae") ) %>%
# dplyr::filter(Period %in% c("SteadyState", "Inhibition", "Recovery/Foaming", "Recovery/Feedingpause",
# "Glycerol"))
ps.summary_sel <- ps.summary %>%
dplyr::filter( str_detect(Family, "Nocard") | str_detect(Family, "Microthrixaceae") | str_detect(Family, "Corynebacter") |
str_detect(Family, "Mycobacter") | str_detect(Family, "Tsukamurellaceae") |
str_detect(Genus, "Tetrasphaera") | str_detect(Genus, "Methanospirillum")
)
library(ampvis2)
psobject <- ps.temp
#set.seed(12345)
#psobject <- phyloseq::rarefy_even_depth(psobject, replace = FALSE)
psobject <- subset_samples(psobject, sample_data(psobject)$week %in% c(10, 11,12, 13,14,15,16))
psobject <- prune_taxa(taxa_sums(psobject) > 0, psobject) #minimum reads per feature
taxtable <- data.frame(tax_table(psobject))
ASVtable <- data.frame(otu_table(psobject))
metaobject <- data.frame(sample_data(psobject)) %>%
rownames_to_column("#OTU ID")
amp <- amp_load(otutable = ASVtable,
metadata = metaobject,
taxonomy = taxtable)
## Warning: Could not find a column named one of "OTU, ASV, #OTU ID" in otutable.
## Using row names as OTU ID's.
## Warning: Could not find a column named one of "OTU, ASV, #OTU ID" in taxonomy.
## Using row names as OTU ID's.
filtervec <- ps.summary_sel$Genus %>% unique()
amp <- amp %>% amp_subset_taxa( normalise = TRUE, tax_vector = filtervec )
## 8539 OTUs have been filtered
## Before: 8777 OTUs
## After: 238 OTUs
p1 <- amp_heatmap(amp, tax_aggregate = "Genus", tax_add = c("Family"),
group_by = "daysop",
tax_show = 30, facet_by = "SludgeType",
showRemainingTaxa = TRUE, normalise = FALSE,
plot_values_size = 3,
min_abundance = 0.001)
p1 <- p1 + theme(strip.text.x = element_text(size = 8, colour = "black", angle = 90),
axis.text.x = element_text(size = 13))
ggsave(plot = p1, "./Figures/Sup_heatmap_filamen_hydrop.png", height=15, width=25, units='cm', dpi=300)
## Warning in scale_fill_gradientn(colours = color.pal, trans = plot_colorscale, :
## log-10 transformation introduced infinite values.
Supp. figure S21 heatmap filaments kraken2
### ampviz heatmap
ps.temp <- pskr
sample_sums(ps.temp)
## S561toS566_AD.AD6 S1001_3_4_AD1_3_4 S1002_5_6_AD2_5_6 S1147_49_50_AD1_3_4
## 2168798 2163187 2756729 2696911
## S1148_51_52_AD2_5_6 S1154_57_58_AD2_5_6 F110 F124
## 2829157 2339807 2756169 2864909
## F179 F5 F55 S119.full
## 2703627 2772044 2692623 1627022
## S1273.ADFull S1554.ADFull S1813.ADFull S370.full
## 1738128 1337619 1489269 1455364
sample_data(ps.temp)$SludgeType <- factor(sample_data(ps.temp)$SludgeType , levels = c("Control","Treatment","Foam","Full-scale", "PSTWAS"))
ps.summary <- ps.temp %>% phyloseq::psmelt() %>%
dplyr::select(OTU, Date, SludgeType, week, Abundance, Kingdom, Phylum, Class, Order, Family, Genus, Species) %>%
dplyr::filter(Abundance > 0)
#ps.summary_sel <- ps.summary %>%
# dplyr::filter( str_detect(Family, "Nocard") | str_detect(Family, "Microtrichaceae") | str_detect(Family, "Corynebacter") |
# str_detect(Family, "Mycobacter") | str_detect(Family, "Tsukamurellaceae") ) %>%
# dplyr::filter(Period %in% c("SteadyState", "Inhibition", "Recovery/Foaming", "Recovery/Feedingpause",
# "Glycerol"))
ps.summary_sel <- ps.summary %>%
dplyr::filter( str_detect(Family, "Nocard") | str_detect(Family, "Microthrixaceae") | str_detect(Family, "Corynebacter") |
str_detect(Family, "Mycobacter") | str_detect(Family, "Tsukamurellaceae") )
library(ampvis2)
psobject <- ps.temp
#set.seed(12345)
#psobject <- phyloseq::rarefy_even_depth(psobject, replace = FALSE)
psobject <- subset_samples(psobject, sample_data(psobject)$week %in% c(10, 11,12, 13,14,15,16))
psobject <- prune_taxa(taxa_sums(psobject) > 0, psobject) #minimum reads per feature
taxtable <- data.frame(tax_table(psobject))
ASVtable <- data.frame(otu_table(psobject))
metaobject <- data.frame(sample_data(psobject)) %>%
rownames_to_column("#OTU ID")
amp <- amp_load(otutable = ASVtable,
metadata = metaobject,
taxonomy = taxtable)
## Warning: Could not find a column named one of "OTU, ASV, #OTU ID" in otutable.
## Using row names as OTU ID's.
## Warning: Could not find a column named one of "OTU, ASV, #OTU ID" in taxonomy.
## Using row names as OTU ID's.
filtervec <- ps.summary_sel$Genus %>% unique()
amp <- amp %>% amp_subset_taxa( normalise = TRUE, tax_vector = filtervec )
## 11684 OTUs have been filtered
## Before: 12380 OTUs
## After: 696 OTUs
p1 <- amp_heatmap(amp, tax_aggregate = "Genus", tax_add = c("Family"),
group_by = "daysop",
tax_show = 30, facet_by = "SludgeType",
showRemainingTaxa = TRUE, normalise = FALSE,
plot_values_size = 3,
min_abundance = 0.001)
p1 <- p1 + theme(strip.text.x = element_text(size = 8, colour = "black", angle = 90),
axis.text.x = element_text(size = 13))
ggsave(plot = p1, "./Figures/Sup_heatmap_filamen_hydrop_kraken.png", height=15, width=20, units='cm', dpi=300)
## Warning in scale_fill_gradientn(colours = color.pal, trans = plot_colorscale, :
## log-10 transformation introduced infinite values.
Supp. figure S22 top 30 species Kraken2/Bracken
psobject <- pskr
psobject <- phyloseq::subset_taxa(psobject, Kingdom %in% c("Archaea", "Bacteria"))
psobject <- prune_taxa(taxa_sums(psobject) > 0, psobject)
set.seed(1000)
psobject <- phyloseq::rarefy_even_depth(psobject)
## You set `rngseed` to FALSE. Make sure you've set & recorded
## the random seed of your session for reproducibility.
## See `?set.seed`
## ...
## 1807OTUs were removed because they are no longer
## present in any sample after random subsampling
## ...
metaobject <- data.frame(sample_data(psobject)) %>%
rownames_to_column("#OTU ID")
#metadata$`#OTU ID` <- str_replace_all(metadata$`#OTU ID`, "-", ".")
taxtable <- data.frame(tax_table(psobject))
ASVtable <- data.frame(otu_table(psobject))
amp <- amp_load(otutable = ASVtable,
metadata = metaobject,
taxonomy = taxtable)
## Warning: Could not find a column named one of "OTU, ASV, #OTU ID" in otutable.
## Using row names as OTU ID's.
## Warning: Could not find a column named one of "OTU, ASV, #OTU ID" in taxonomy.
## Using row names as OTU ID's.
p1 <- amp_heatmap(amp, tax_aggregate = "Species", tax_add = c("Phylum"), tax_show = 30, showRemainingTaxa = TRUE, normalise = FALSE, group_by = "SludgeType2")
p2 <- amp_heatmap(amp, tax_aggregate = "Species", tax_add = c("Phylum"), tax_show = 30, showRemainingTaxa = TRUE, normalise = TRUE,
group_by = "daysop", facet_by = "SludgeType") +
theme(
axis.text.x = element_text(size = 13)
)
#g2 <- ggarrange(p1, p2)
p2
## Warning in scale_fill_gradientn(colours = color.pal, trans = plot_colorscale, :
## log-10 transformation introduced infinite values.
# sum(taxa_sums(psobject)) #15,219,740 reads
# 4,220,112 reads remaining
# sum(taxa_sums(psobject)) - 4220112 # 10999628 reads represented in figure
# 10999628 / 15219740 # 72.3%
# sum(taxa_sums(psobject)) #9,930,576 23 July 24, including PSTWAS, full-scale, rarefied
# 1982212 remaining
# sum(taxa_sums(psobject)) - 1982212 #7948364 reads represented in figure
# 7948364 / sum(taxa_sums(psobject)) #80%
ggsave(plot = p2, "./Figures/Sup_heatmap_kraken_RMITADs.png", height=15, width=31, units='cm', dpi=300)
## Warning in scale_fill_gradientn(colours = color.pal, trans = plot_colorscale, :
## log-10 transformation introduced infinite values.
Supp. figure S23 top 30 species
ps.onp <- ps.onp.ncbi.kraken
#summary(sample_sums(ps.onp))
set.seed(1000)
psobject <- phyloseq::rarefy_even_depth(ps.onp, replace = TRUE)
## You set `rngseed` to FALSE. Make sure you've set & recorded
## the random seed of your session for reproducibility.
## See `?set.seed`
## ...
## 2041OTUs were removed because they are no longer
## present in any sample after random subsampling
## ...
# ps.onp.ncbi (minimap) phyloseq-class experiment-level object
# otu_table() OTU Table: [ 4866 taxa and 8 samples ]
# sample_data() Sample Data: [ 8 samples by 19 sample variables ]
# tax_table() Taxonomy Table: [ 4866 taxa by 8 taxonomic ranks ]
#ps.onp.ncbi.kraken - phyloseq-class experiment-level object
#otu_table() OTU Table: [ 5166 taxa and 8 samples ]
#sample_data() Sample Data: [ 8 samples by 19 sample variables ]
#tax_table() Taxonomy Table: [ 5166 taxa by 8 taxonomic ranks ]
#ps.onp.st8 (kraken) - phyloseq-class experiment-level object
#otu_table() OTU Table: [ 2472 taxa and 8 samples ]
#sample_data() Sample Data: [ 8 samples by 19 sample variables ]
#tax_table() Taxonomy Table: [ 2472 taxa by 8 taxonomic ranks ]
taxtable <- data.frame(tax_table(psobject))
ASVtable <- data.frame(otu_table(psobject))
metaobject <- data.frame(sample_data(psobject)) %>%
rownames_to_column("#OTU ID")
metaobject$Sample_name <- str_replace_all(metaobject$Sample_name, c("CSTR1"="C1",
"CSTR2"="T1",
"CSTR3"="C2",
"CSTR4"="C3",
"CSTR5"="T2",
"CSTR6"="T3") )
amp <- amp_load(otutable = ASVtable,
metadata = metaobject,
taxonomy = taxtable)
## Warning: Could not find a column named one of "OTU, ASV, #OTU ID" in otutable.
## Using row names as OTU ID's.
## Warning: Could not find a column named one of "OTU, ASV, #OTU ID" in taxonomy.
## Using row names as OTU ID's.
#p1 <- amp_heatmap(amp, tax_aggregate = "Genus", group_by = "SludgeType", tax_show = 30, showRemainingTaxa = TRUE, normalise = FALSE)
p <- amp_heatmap(amp, tax_aggregate = "Species",
tax_add = c("Phylum"),
facet_by = "SampleType",
group_by = "Sample_name",
plot_values_size = 3,
tax_show = 50, showRemainingTaxa = TRUE, normalise = TRUE)
#ggsave(polot = p, ggsave("./Figures/heatmap_16SFL.png", height=22, width=22, units='cm', dpi=300))
Supp. figure EC (not included)
firstdate <- 47 # selected periods for general manuscript
lastdate <- 124
weeky <- c(-0.5)
labely <- c(0.5)
segmenty <- c(0)
generate_vector <- function(x, y) {
result <- c()
while (x <= y) {
result <- c(result, x)
x <- x + 7
}
return(result)
}
generate_vectorday <- function(x, y) {
result <- c()
while (x <= y) {
result <- c(result, x)
x <- x + 1
}
return(result)
}
foamdata %>%
dplyr::filter(SludgeType %notin% c("Full-scale")) %>% # filter columns out
ggline(x = "DayofOp",
y = "EC",
color = "SludgeType",
legend = "right",
add = "mean_se",
add.params = list(width = 1)) +
# scale_color_manual(values = cols) + # create a cols vector for custom colours
theme_light() +
scale_color_manual("Sludge Supernatent", values = cols) +
theme(
axis.text.x = element_text(angle = 0, hjust=0.5),
legend.position = c(0.90, 0.825),
legend.background=element_rect(fill = alpha("white", 0.5)),
legend.ticks.length = element_blank(),
legend.title = element_blank() ) +
scale_x_continuous(breaks = generate_vector(firstdate,lastdate),
minor_breaks = generate_vectorday(firstdate,lastdate),
limits = c(firstdate, lastdate) ) +
annotate("text",x = 51, y = weeky, label = "W08 ", size = 2.5) + #Week 08
annotate("text",x = 58, y = weeky, label = "W09 ", size = 2.5) + #Week 09
annotate("text",x = 65, y = weeky, label = "W10 ", size = 2.5) + #Week 10
annotate("text",x = 72, y = weeky, label = "W11 ", size = 2.5) + #Week 11
annotate("text",x = 79, y = weeky, label = "W12 ", size = 2.5) + #Week 12
annotate("text",x = 86, y = weeky, label = "W13 ", size = 2.5) + #Week 13
annotate("text",x = 93, y = weeky, label = "W14 ", size = 2.5) + #Week 14
annotate("text",x = 100, y = weeky, label = "W15 ", size = 2.5) + #Week 15
annotate("text",x = 107, y = weeky, label = "W16 ", size = 2.5) + #Week 16
annotate("text",x = 114, y = weeky, label = "W17 ", size = 2.5) + #Week 17
annotate("text",x = 121, y = weeky, label = "W18 ", size = 2.5) + #Week 18
annotate("text",x = 54, y = labely, label = "Steady state", size = 2.5) +
annotate("segment", y = segmenty, yend = segmenty, x = 47, xend = 62,
arrow = arrow(ends = "last", angle = 90, length = unit(.1,"cm"))) +
# annotate("text",x = 65, y = labely, label = "Glycerol ", size = 2.5) +
annotate("segment", y = segmenty, yend = segmenty, x = 62, xend = 67,
arrow = arrow(ends = "both", angle = 90, length = unit(.1,"cm"))) +
annotate("text",x = 74, y = labely, label = " Inhibition", size = 2.5) +
annotate("segment", y = segmenty, yend = segmenty, x = 67, xend = 82,
arrow = arrow(ends = "both", angle = 90, length = unit(.1,"cm"))) +
annotate("text",x = 86, y = labely, label = "Recovery", size = 2.5) +
annotate("segment", y = segmenty, yend = segmenty, x = 82, xend = 90,
arrow = arrow(ends = "both", angle = 90, length = unit(.1,"cm"))) +
#annotate("text",x = 94, y = labely, label = "Foaming ", size = 2.5) +
annotate("segment", y = segmenty, yend = segmenty, x = 90, xend = 98,
arrow = arrow(ends = "both", angle = 90, length = unit(.1,"cm"))) +
annotate("text",x = 101, y = labely, label = " Recovery", size = 2.5) +
annotate("segment", y = segmenty, yend = segmenty, x = 98, xend = 106,
arrow = arrow(ends = "both", angle = 90, length = unit(.1,"cm"))) +
annotate("text",x = 115, y = labely, label = "Steady state post foam", size = 2.5) +
annotate("segment", y = segmenty, yend = segmenty, x = 98, xend = 124,
arrow = arrow(ends = "first", angle = 90, length = unit(.1,"cm"))) +
annotate("text",x = 85, y = 4, label = "Feeding paused", size = 4, angle = 90, alpha = 0.5) +
annotate("rect", ymin = 0, ymax = 8, xmin = 82, xmax = 89, fill = "blue", alpha = .1) +
annotate("text",x = 64, y = 4, label = "Glycerol", size = 4, angle = 90, , alpha = 0.5) +
annotate("rect",ymin =0, ymax = 8, xmin = 62, xmax = 67,
fill = "blue", alpha = .1) +
annotate("text",x = 93.5, y = 4, label = "Foam observed", size = 3.5, angle = 90, , alpha = 0.5) +
annotate("rect", ymin = 0, ymax = 8, xmin = 91, xmax = 96,
fill = "blue", alpha = .1) +
ylab("EC (S/cm)") +
xlab("Day of operation")
ggsave("./Figures/Sup_lineplot_ec_daysop.png", height=2.5, width=6, units='in', dpi=300)